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ABSTRACT 



This thesis presents a model which simulates the scattering from a fluid 
loaded I-beam and the resultant behavior due to fluid-structure interaction. 
Chapter I gives an overview of the problem and describes the characteristics 
of the solid and fluid, the aspects of periodicity, boundary conditions and the 
coupling of the two media. 

The governing equations of motion are scaled in Chapter II. In Chapter 
III, the finite-difference formulae for these equations are derived, as is the 
non-local radiation boundary condition. Difference formulas for typical 
boundary points of the solid and corner nodes are also derived. All finite 
difference formulae used are presented in Appendix C. Chapter IV contains 
numerical results. Conclusions are drawn and areas of the problem that 
would require further study are in Chapter V. 
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I. PROBLEM OVERVIEW AND DOMAIN DESCRIPTION 



The problem we consider is where an acoustic pressure wave emitted by 
an active sonar impinges on a submarine. In our case this is a double-hulled 
submarine. Active sonar works on the principle of detecting the reflection off 
a solid object of an acoustic pressure wave emitted by a source, by which one 
can calculate distance and direction to that object as in Figure 1, but instead of 
concentrating on the reflected wave we look at the interaction between the 
structure and the incoming wave. This generates scattered pressure waves. 
The scattered pressure waves include waves which decay as they travel 
through the fluid (evanescent modes), and waves which do not (propagating 
modes). Since propagating modes do not decay they can be detected. The 
main thrust of this thesis is to determine the characteristics of the propagating 
modes, such as amplitude and energy. The optimum situation would be to 
perturb the double-hulled structure at a resonance frequency which would 
increase the amplitude of the propagating modes making them easier to 
detect. We investigate the steady state behavior of the propagating modes and 
the resultant shear strain field in the I-beam. At steady state the 
characteristics of the propagating modes stabilize which might be used as an 
acoustic signature of the structure. In addition when the shear strain field 
reaches steady state its value throughout the solid may be used to isolate areas 
of the I-beam where large stresses occur. This information could be useful in 
the design of double-hulled vessels. 
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Figure 1. Physical Situation 
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To reduce the model to one where we can apply numerical techniques to 
simulate its behavior we must make simplifying assumptions which are 
discussed below. 



A. ASSUMPTIONS 

• The source of the pressure waves is far enough away so that the 
impinging wave can be approximated by a normally incident plane 
wave. 

• There are no other sources present such as might occur with bottom 
bounce, surface reflection and the like. 

• The area of interest — AOI as shown in Figure 2 where the wave 
impinges is where the inner and outer hulls are joined or connected by 
a supporting spar forming an I-beam shaped domain. 

• The I beam shaped domain is a uniformly continuous linear isotropic 
elastic medium with no cracks, welds or other deformities. 

• The dimensions of our AOI are such that any curvature of the surfaces 
can be ignored. 

• The center spar or support beam occurs at regular intervals through the 
structure as shown in Figure 3 allowing us to truncate the domain to 
the left and right of center spar using periodic boundary conditions. 

• The incident wave does not displace the submarine. 

• The cavities A and B in Figure 2 enclosed by the inner and outer hulls, 
and the center spar is void and contain no sources. 

• We are only interested in displacement in the x and y directions as 
given in Figure 2 which reduces our problem to two dimensions. 

• The fluid is seawater, the solid is steel. 

Our model is now reduced to one where we have two coupled media — fluid 
and solid, and the characteristics of each will now be discussed in turn. 
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Figure 2. Area of Investigation 
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A. THE SOLID 

A cross-section of the double hull of length 8a is shown in Figure 3. Note 
a is a scaling constant for distance. The domain is essentially infinite in 
extent in the x (length) and z (depth) directions. Since it is of uniform shape 
in the z direction, we can neglect this coordinate and reduce our problem to 
two dimensions. Our area of investigation is outlined in Figure 3, and we are 
now faced with the problem of providing boundary conditions in the x 
domain. 




Figure 3. Cross-section of Submarine Hull 



We have a normally incident plane wave impinging over one boundary 
of the domain. Thus the pressure is the same at all points (there is no phase 
shift) along the fluid /structure boundary. Due to the periodic nature of the 
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structure we expect the solid to behave the same over given intervals of 2a, 
that is 

f(x,y,t) = f(x + 2aN,y / t) (1) 



where N is an integer and / a function describing the state of the solid such as 
displacement, velocity and so on. Using this periodicity we can truncate our 
domain at a and -a and use the following periodic boundary conditions 

f{-a,y,t) = f{a,y,t) (2) 



and 



f k {-a,y,t) = f k (a,y,t) 



(3) 



where k is a positive integer and can signify the derivative with respect to x, y 
or f, depending on the function f. 

The solid has now been reduced to a two-dimensional linear isotropic 
elastic medium whose governing equation of motion for points interior to 
the domain is given by the plane strain elastic wave equation which is 



^ 3 2 t< d^u 

a? + a/ 



+(A +/i) 



^d 2 u 



aV 



dx 2 dxdy 



= Ps 



d 2 U 

dt 2 



(4) 



^ d 2 V 



d 2 V^ 



dx 2 3/ 



+ (A +/x) 



f d 2 V 



aV 



dx 2 dxdy 



~ Ps 



d^V 

dt 2 



(5) 



u and v are displacements in the lateral (x) and transverse (y) direction, [X and 
A are Lame constants, and p s is the density of the solid. The I-beam is 
deformed due to the imposed fluid pressure. Balancing normal forces at the 
fluid /solid interface we obtain the boundary condition 
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( 6 ) 



Tyy - A 



du dv^\ . dv total 

^: + T- +2 vjc = -P 

\dx oy J ay 



where p total is the total pressure, and Tyy the normal stress component. All 
other boundary conditions are either traction free or periodic and are 
summarized below. 



T xx ~ 0 
T xy ~ ^ 



► on surface DH and El in Figure 4, 



(7) 



%=° 
T xy = 0 



► on CD, EF, GH, IJ, and KL in Figure 4, 



(8) 



T w = -P 

^ xy = ® 



total ' 



on AB in Figure 4, and 



(9) 



periodic 

conditions 



elsewhere- 



( 10 ) 



t XX/ Tyy and r X y are the normal and shear components of stress and are 
represented by 



Tyy - A 



du dv ^ 
dx + dy. 



+ 2ji 



r dv' 



\ d yj 



( 11 ) 



T'XX ~ ^ 



du dv ^ 
dx dy 



+ 2 [i 



%' and 



( 12 ) 



*xy = /* 



du dv ^ 
dy + dx/ 



(13) 
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Figure 4 is the unsealed domain, in which s, and k are constants used to vary 
the size of the cavities A and B of Figure 2. Note: 0 < s < 1, 0<k< 1. 



U 



FLUID 



Al 



I B 




x = -a 



x = 0 



x = a 



Figure 4. Unsealed Domain 



B. THE FLUID 

We are interested in the scattered pressure waves generated by 
perturbations at the fluid/solid interface. The partial differential equation 
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governing the propagation of these waves in the fluid is the linear two- 
dimensional wave equation given by 



1 d 2 p _ d 2 p d 2 p 
c 2 dt 2 dx 2 dy 2 



( 14 ) 



where cj is the speed of sound in the fluid. 

These perturbations arise when a normal plane wave of magnitude P and 
frequency co impinges on the surface of the solid, causing it to deform. Since 
the pressure wave is of normal incidence (there is no phase shift at points 
along the surface of the solid) and the solid is a periodic structure we can 
expect the behavior of the scattered pressure waves to be the same in given 
intervals of 2a, thus we can use periodic boundary conditions in the same 
manner as was used for the solid. 

Ideally, if the solid were acoustically hard our total pressure would only 
have two components, incident and reflected. In reality the solid is perturbed 
by the incident wave, generating scattered pressure waves which propagate 
out into the fluid domain. The total pressure in the fluid can be represented 

ty 

ptotal = ^incident + ^reflected + ^scattered (14) 

where ^incident = e" %k P~ Xi ) ^reflected = e lk fV~ lt and ^scattered i s t o be 
determined. To avoid cavitation at the fluid-solid boundary we employ the 
inviscid form of the Navier-Stokes equation which we refer to as the 
compatibility condition and is given by 
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( 15 ) 



where pf is the density of the fluid. Note that only the scattered term of the 
pressure is included in this equation since 



The scattered pressure waves generated at the fluid/solid interface are 
composed of propagating (non-decaying) and evanescent (decaying) modes 
which must be allowed to propagate off to infinity in the positive y direction. 
To do this we must employ a non-local radiation boundary condition whose 
implementation will be discussed in greater detail in the discretization 
section of this paper. 




( 16 ) 



thus 




( 17 ) 
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n. NON-DIMENSIONALIZATION OF THE GOVERNING EQUATIONS 
OF MOTION AND BOUNDARY CONDITIONS 



The problem addressed models the propagation of a scattered pressure 
wave in two dimensions. This is described by the two-dimensional scalar 
wave equation 



1 d 2 p = d 2 p d 2 p 
cj dt 2 dx 2 + dy 2 



(18) 



where p is the pressure and Cf the speed of sound in the fluid. To facilitate 
implementation and to free ourselves from the requirement of using a given 
system of measurement such as metric or imperial we scale or non- 
dimensionalize Equation 18 as follows. Let 



1 d 2 p _ d 2 p d 2 p 
c 2 dt 2 dx 2 d y 2 



(19) 



represent the unsealed wave equation. We now use the following 
relationships: 



t = cot; 





( 20 ) 



Here co is the scaling constant for time and the frequency of the incident plane 
wave, a is the half length of the I-beam and the scaling constant for distance 
and P is the scaling constant for pressure. Note that when taking derivatives 
with respect to x we get 
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and 



( 21 ) 



d__ \d_ 
dx a dx 



dx 2 a 2 dx 2 



This also holds true for derivatives with respect to y, and a similar relation 
results when taking derivatives with respect to t. With the above 
relationships, Equation 19 is now written as 



0) 2 P d 2 p _ P d 2 p P d 2 p 

cj 3f 2 " a 2 3x 2 + a 2 dy 2 ' 



( 22 ) 



Cancelling common factors and multiplying Equation 22 by a 2 reduces it to 



co 2 fl 2 cPp _ 5 2 p d^p 
c 2 dt 2 dx 2 dy 2 



(23) 



Defining 




(24) 



Equation 23 is now written as 



k 2^P = *y + ^P 



''a? 



dx z dy 1 



(25) 



The elastic wave equation, which is a vector wave equation is given by 

f ^2 "\2 

d u 3 m 



dx 2 dy 2 



+ (A +p) 



( -\2 <2 \ 

d U d V 

+ 



dx 1 dxdy 



d 2 u 

- ps d t 2 



(26) 



^ d^v_ d 2 v ' 
dx 2 dy' 



+ (A+p) 



( o <2 \ 

d V dl 
dy 2 dxdy 



d 2 V 

~ Ps tf' 



(27) 
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The following relationships allow us to write Equations 26 and 27 in a more 
convenient form. 



— = Cj and ^- = c 2 L . 
Ps Ps 



(28) 



The constants ci and cj are the lateral and transverse velocities of the solid. 
The unsealed forms of Equations 26 and 27 are now written as 

d 2 u d 2 u ) /, 7 \( d 2 u ^ ^ 



4 



dx 2 df 



+ 



(4~4) 



d 2 v ^ 



dx 2 



dxdy 



d z u 



dt 



T 2 



(29) 



4 



-.2— ■'v2 — '\ 

d V d V 



dx 2 dy 



—2 



( 4-4 ) 



d 2 V + d 2 u ^ 
v 9y 2 dxdy j 



d 2 V 



dt 
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(30) 



We use the same scaling relationships as before for x, y and t (Equation 20) as 
well as 



u v 

u = ~p and v =jj , 



(31) 



which are the scaling relationships for the displacement, and rewrite 
Equations 29 and 30 as 



4 



D d 2 u 
a 2 dx 2 



D d 2 u 

Idy 2 



( 4-4 ) 



D d 2 u D 



a 2 dx 2 



+ ■ 



d 2 v 



\ 



dxdy 



= Do 



2 d 2 U 

4? 



(32) 



4 



D d 2 v D d 2 v 



dx' 



a 2 dy 2 



(4-4) 



D d 2 v D d 2 u 



a 2 dy 2 



+ ■ 



, 2 dxdy 



n 2 ° v 
= Deo — T . 

dt 2 



(33) 



Defining 




and 




(34) 
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and dividing Equations 32 and 33 through by co 2 and cancelling common 
factors they reduce to 



, /O ^2 \ f i i Y ->2 -\2 \ -\2 

1 d u d u 1 1 d u d v d u 



k T 



dx 2 dy z 



\ k l k tJ 



dx 2 dxdy 



dt* 



(35) 



ki 



r d 2 v + d 2 v ' 

a? dp 



j T_ 



(<2 



2 .. ^ 



+ 



ay 2 a*ay 



a^£ 

at 2 ' 



(36) 



Collecting like terms gives us the final form 



1 d 2 u 1 d 2 u 
k T dy 2 



k } dx 2 



V > 



l2 l 2 

K k l k t j 



dy_ 

dxdy 



d^u 
a t 2 



(37) 



1 d 2 v 1 d 2 v 

+ —fT 7T + 



kl dx 2 



2 -v,2 



k L dy 



V ->2,. ^ 



V*L 



ju2 

k Tj 



d^U 

dxdy 



d 2 v 

a?' 



(38) 



The surfaces in contact with free space are stress free. (The fluid/solid 
boundary is dealt with separately). Therefore the stresses t xx and t X y on 
surfaces El and DH and t yy and x X y on surfaces CD, EF, GH, IJ and KL of 
Figure 5 are zero. 

These components of the stress tensor can be written 



X xy ~ A* 



du + dv^ 
dy dx 



= 0 



(39) 



T xx =l 



du dvP 
\dx + dy 



_ du n 

+ 2,- = 0 



(40) 



%- A 



k du dv^ 



dx + dy 



_ dv . 
+ 2^1 — = 0 

dy 



(41) 
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Figure 5. Scaled Domain 



where x, y , u and v are the unsealed components of distance and 
displacement and fj. and X are the Lame constants. We will now scale 
Equations 39 through 41 in turn. 

For Equation 39 using the same scaling relationships as before (see 
Equations 20 and 31) gives 
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(42) 






D du + D dv 
a dy a dx 



= 0 . 



Cancelling common factors reduces it to 



du dv 

T *y = dy + ~d£ = °- 



(43) 



For Equation 40 we first divide through by p s the density of the solid, and 
using the same scaling relationships we get 



A ( D du Ddv 



Ps 



a dx a dy ) p s a dx 



2 p D du _ 
+ — — — = 0 . 



(44) 



We know that 



p/ps = c T 



(45) 



Similarly it can be shown that 



A 2 2 

p s - c l- *T- 



(46) 



Using Equations 45 and 46, substituting into Equation 44, and cancelling 
common factors gives 



(^- 2 4 | + | 



dx 



(47) 



or 



r r 2 

C L 



-2 






c T 



dx dy 



„du _ 
+ 2 — = 0 . 

dx 



( 48 ) 



Using the previous definitions of ki and kj it can be shown that 
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( 49 ) 



±J± 

c\ kl ' 



and when used in Equation 48 it reduces to our final form 



$du 

kite 



(*l ) 



= o. 



(50) 



Following the same procedure for Equation 41 its final form is 



K*l 



dx kl 



(51) 



At the interface the fluid cannot exert a force tangential to the boundary, 
hence the shear component of stress is zero there, and is given by 

( du dv ^ 



*xy=H 



__ + _ 

dx dy 



= 0- 



(52) 



For the normal component of stress at the interface we refer back to Equation 
6, Section B and Equation 15 in Section C of Chapter I to see that the normal 
component of stress is given by 



x yy 



= X 



'du dv N 
^ dx dy j 



+ 2 ^ 



du 

dx 




+ p R + p S ), 



(53) 



where the superscripts I, R and S signify incident, reflected and scattered 
pressures. The incident and reflected pressures are given by 

p z =e -»/y-« and p K =e */y-« 

x xy is scaled the same as before. To scale x yy properly we will need the 
compatibility condition, whose unsealed form is 
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( 54 ) 



— s 



dp 






d 2 v 

y = o Pf di 2 ' 



We refer back to Equation 17, Chapter I, Section C to see that only the 
scattered pressure need be considered here. 

Scaling Equation 54 we get 



P dp s 
a dy 



y = o 





(55) 



or 



-P 




= 0 



■PfCQ 



-aD 



d 2 V 



dt 



2 ' 



(56) 



Equation 56 gives us a convenient choice for the scaling constant for pressure 
of 



P = (0 2 aDpj 



(57) 



and when substituted back into Equation 56 cancels as a common factor to 
give 



dp 5 



a y 



y = 0 



d 2 V 

a?' 



(58) 



Dividing Equation 53 by p s and scaling gives 



W 

Ps a 



du | dv 
dx dy 



2 uD dv -P t i 

+ -P—— = — Ip 1 +p 
p s a dy p s \ 



R 




Substituting the value of P from Equation 57 we obtain 



( 59 ) 
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( 60 ) 



AD 

Ps a 



du dv 
dx dy 



2 uD dv Pf 2 / / r c\ 

+ — = — -coaDip +p +v } 

Ps a Ps ^ ' 



y = o 



or 



0£_ 2 



,2/ / r s\ 

— + — +2 — = -ekj\p +p K + tr 
dx dy) dy T \ F y F > 



A 



y = o 



(61) 



where e = pf/ps- 

Evaluating the incident and reflected pressures at y = 0 we get 
pi = pR = e~ il , and when substituted into Equation 61 gives 



fa 
%-2 
k } 



^fdu dv 

— + 



\ K L 



dx dy 



+ 2— = -lekje lf - ekjp s . 

dy 



(62) 
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III. DERIVATION OF THE FINITE DIFFERENCE FORMULAE FOR THE 

GOVERNING EQUATIONS 

Throughout the derivation of the finite difference approximations we 
identify gridpoints of the scaled domain with the subscripts i and j and 
superscript n. Variation in the x direction is denoted by i, 1 < i < N, where N 
is the total number of subdivisions (since we are using a square grid the 
number of subdivisions in the y direction is the same), variation in y with j, 1 
< j <N, and time with n, 0 < n < <». Lower case indices denote varying 
quantities, while upper case letters denote fixed quantities. would be the 
value of the scattered pressure for all values of i at the grid level y = KAy, thus 
p” K is a vector, while is an element. 

As was discussed above we use the letter i to denote variation in the x 
direction. This is a common convention and we do not wish to deviate from 
it. To avoid confusion with the complex quantity, /-I , also commonly 
denoted by i, we state the following rule, that whenever the letter i appears as 
a superscript it denotes the complex quantity yj-l, and when appearing as a 
subscript it is an index denoting variation in the x direction. 

A. FINITE DIFFERENCE APPROXIMATIONS FOR THE EQUATIONS 

GOVERNING THE BEHAVIOR OF THE FLUID 

The scaled domain of consideration as shown in Figure 6 has periodic 
boundary conditions applied at x = 1 and -1 and a radiation boundary 
condition for the propagating modes at y = 2. To model the two-dimensional 
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wave equation numerically we must derive an equivalent finite difference 
formula, from which we can solve for the pressure for all subsequent time 
levels as follows: 



2 a 2 ?? _ a 2 ? a 2 p 

f dt 2 ~dx 2 dy 2 



(63) 



is the two-dimensional wave equation as derived in Chapter II Section A. 
Using central difference approximation for all second order partial 
derivatives the equivalent finite difference equation for Equation 63 is 

At 2 




(64) 



where h = Ax = Ay is the step size and At is the increment in time. The 

truncation error for Equation 64 is 0(/i 2 ) in space and 0(At 2 ) in time. Solving 
for p” ; +1 explicitly we have 




Letting p 2 



Equation 65 can be written as 

kjh 2 




( 66 ) 



The Von Neumann stability criterion 

1 




( 67 ) 
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must be satisfied to ensure the stability of Equation 66. 1 

Special attention must be paid when applying the wave equation along 
the boundaries at x = 1 and -1 for it is here that we make use of the periodic 
boundary conditions. Applying Equation 66 at x = -1, i = 0 and y = jh (see 
Figure 7) yields 




This requires the value of p at the point (-1, jh, nAt ) which lies outside the 
domain. By using the periodic boundary conditions 

p( I, y, t ) = pi-1, y, t ) (69) 



dp 

dx 



= dp 

i+W) dx 



(-W)' 



(70) 



we can substitute the value of ((N-l)h, jh, nAt ) (where N is the total number 
of subdivisions in the x direction) for the value at (-1 ,jh, nAt), allowing us to 
evaluate the wave equation at the boundary. 



B. APPLICATION OF THE RADIATION BOUNDARY CONDITION IN THE 
NUMERICAL SCHEME 

As was mentioned at the end of Chapter I (Section C), we apply a non- 
local radiation boundary condition (referred to as nlrb ) to the fluid to 
simulate an infinite domain in the positive y direction. Our domain is 
truncated at x=l and -1, forcing our fluid domain to act as a waveguide. The 
scattered pressure can be represented as a series of plane waves which take the 
form 



J For treatment of the von Neumann Stability Criterion see Appendix A. 
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Figure 7. Left Boundary for the Fluid 
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e '{rkx+Pky - 0 



where 




for propagating modes 
for evanescent modes 



(71) 



and % = kn. We note that when fa = i^jYk -kj Equation 71 yields 

e -Pky+Kr kx-t) 

which is an exponentially decaying quantity for positive values of y. This fact 
allows us to neglect evanescent modes when applying the radiation boundary 
condition, at y = 2. 

The total scattered pressure is given by 

(72) 

k = —°° 



This is composed of propagating and evanescent modes. Far from the fluid 
solid interface where only the propagating modes are assumed present (for 
reasons given above) the scattered pressure is written 

p= (73) 

k=-M 

where M is the total number of modes (positive or negative) under 
consideration. At the boundary y = 2, we apply the nlrb operator (Scandrett 
and Kriegsmann, 1992, unpublished paper). 

(74) 

to the individual modes of the scattered pressure of Equation 73 which yields 
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M 



k=-M 



M 



BCp)- I £ ft(|(V (rH+A! '-' ) 



*=-M 



(75) 



Equation 75 reduces to 



M 



B(p)= S(iA-/AV (r ‘ JI+fey - ,) . 

k=-M 



(76) 



The right-hand side of Equation 76 is identically zero. The boundary operator 
has annihilated the propagating modes and since the evanescent modes are 
assumed to have negligible magnitude there, any scattered pressure waves 
reaching the boundary experience no reflection, simulating an infinite 
domain in the positive y direction. 

To apply the nlrb operator at the boundary y = 2, we rewrite Equation 75 as 



dp 



M 



y=Jh 

t-T 



1 hjl« k (TY (Mh+nx) Jso 



k=-M 



(77) 



where we evaluate the pressure at a constant time T and along the boundary 
y = ]h (i.e. at y = 2). We incorporate ajt and e~ iT as ake~ iT and define this to be a 



new constant ajt(T). Note that ||«A:(^)|| = || fl jt|| since 




= 1 for all values of T. 



ock(T) is unknown so we must derive an alternative expression to be able to 
evaluate Equation 77. The k th propagating mode can be written as 

v,=aJ^P k y-t) (78) 



and when evaluated at y=Jh and at t = T and employing our new constant 
ock(T ) we obtain 
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Pk y=]h = a k {T)e'M h c' y k x . 
t=T 



( 79 ) 



To isolate a k (T)e^ k ^ we multiply Equation 79 by e * ykX (we apply 



orthogonality) and integrate over the x domain to obtain 



l 



jp(x,Jh,T)e iykX dx = 2aj c (T)e t ^ k ^ h 



(80) 



-1 



or 




(81) 



where % is a dummy variable of integration. We substitute this value of 



which allows us to apply the nlrb at the boundary y = 2 and from Equation 81 
we will be able to evaluate the amplitudes of the propagating modes of the 
scattered pressure. 

We now proceed with deriving the finite difference approximation for 
the radiation boundary condition. Using a central difference approximation 



a k (T)e^ k ^ back into Equation 77 to obtain 




(82) 




t=T 



, Equation 82 is written as 



( 83 ) 



n n + n n 

Kj + 1 + Pi,/-1 
2h 




The trapezoidal rule for integration is 

\ a f( x Y x - ^ (/l + Vl + 2 /3 • • • +2 /n-l + fn ) t 84 ) 

and is used for the integral in Equation 83, however it can be more compactly 
expressed as 



fl 



J ' b 'f(x)dx = h'Y J 8 r f r where 8 r = 



r= 1 



— r = 1 or I 
2 



1 elsewhere. 



(85) 



where 1 is the total number of nodes in the x direction. When substituted 

dp 

into Equation 83, using a central difference approximation for (t,, Jh, t ), 
Equation 83 can be written 



,n 






2h 



2 At 



4 

A k=-M 1 r = 1 



Multiplying by , we have 



2 At{„„ 
h : 



I M 



(p"/ + + I IAV 

r=l fc=-M 






( 86 ) 



(87) 



Define A to be a matrix whose i,r entry is given by 

*(■>)= 

k=-M 



( 88 ) 



Upon substitution into Equation 87 we obtain 
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(89) 



2 Attn 



fow-plnbW'-'V j'-o- 



n n n + 1 n- 1 

We note that Pi j+\> Pi j-i >P r j and p r j are all vectors due to the J 



subscript as described at the beginning of this chapter. The radiation 
boundary condition and the wave equation must both be satisfied at the 
artificial boundary y = 2. Applying either of the two conditions will require 
the use of pseudonodes which lie outside the domain. Through the 
combination of the two equations we will be able to eliminate this 
requirement. Reproducing the two-dimensional wave equation and the 
radiation boundary condition, (substitute k for all x indices in the wave 
equation and radiation boundary condition, since these are dummy indices), 
we have 




2At 

h 2 




/+l Pk.J-l 



)* 



A p n k!i 



M j 



= 0 . 



(91) 



Note that both equations are being evaluated at y = 2 and that the vectors 
p %_. j j and p£ +1 j have a circular shift and are of the same dimension as the 

rest of the vector in Equation 91. Collecting like terms in Equations 90 and 91 
we obtain 



At „n At ( t "Ti - 1 

~j^Pk,]+\ -j^Pk,J-\ + Pk,] + Pk,J 



7' 



Ar 

klh 2 



Pk-hJ 




(92) 
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( 93 ) 



l^L„ n — n n + An n+ ^ - An n ~ l - 0 

1.2 Pk,j+\ .2 Pk,J-\ +A Pk,J A Pk,J “ U - 



At 



Adding — j times Equation 93 to Equation 92 yields 
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(94) 



n+1 



Solving for pj. j in Equation 94 we obtain 



n+l 

p k,l 





- 1 - 
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1 + 0 A 




2kj 
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2At 2 n At 2 n 
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kjh 



ph 



At 2 n 
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J--4U 
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n- 1 

Pk,l 
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(95) 



n n 



The terms differenced in x, that is p k _ | j , p^j, Pk+\ j can be more compactly 

expressed as Tp” r where T is a tridiagonal matrix with 2 - — 5 — y on the main 

At 2 k f h 

diagonal and ■ 2 2 on the sub and super diagonals. We must also allow for 

the periodic boundary conditions when constructing T. To do this we replace 



the ( N , 2) and the (1, N- 1) elements of T with „ _ where T is an N by N 

kjh 2 

matrix. The general form of T can be seen from Equation 16 Appendix C. 
Equation 95 is now written as 
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(96) 
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-1“ 


, At 
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2k f J 





2 At 



2 U 2 Pk,]-\ +T Pk,] 



i-*La 



kfh 



\ 



2k 



f 



n n ~ l 

Pk,j 



J 



which satisfies both the wave equation and the radiation condition at the 
boundary. 



C FINITE DIFFERENCE APPROXIMATION FOR THE ELASTIC WAVE 

EQUATION 

To investigate the propagation of disturbances in an / 'I-beam ,/ -shaped 
domain as shown in Figure 5, it will be necessary to apply the elastic wave 
equation to points in the interior of the domain. (The boundaries will be 
dealt with in a separate section.) By only considering displacements in the x 
(lateral) and y (transverse) directions, the problem becomes one of plane 
strain in two dimensions. Reproducing the scaled equations derived earlier 
for motion in the lateral and transverse directions 



1 d 2 u 
k 2 L dx 2 + 



1 d 2 u fj_ 

k$dy 2+ {k 2 L 



kr ) 



d 2 v 



d xdy 



d 2 u 

dt 2 



(97) 



1 d 2 V 1 d 2 V 1 

/4d?' + k[d^ 2 (k! 



i ) a 2 » 
k 2 )dxdy 




(98) 



we use central difference formulas for the partial derivatives in Equation 97 
to get the equivalent finite difference equation 



h 2 k 2 L 



u n _2 U ” +u? 

u i+\,j zw i,/ + u i-l,;J 



klh 2 



, n -2u n +u n 

*i.j+l ZU i,j + 
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K^L k 2 j 



4h 2 [ V M,i + 1 r M,/+ 1 + D i-1,/-1 






p,7 1,7 h ] 



( 99 ) 



The truncation error for Equation 99 is 0(/i 2 ) in space and 0(At 2 ) in time. 



n+l 



Solving for u • j explicitly in Equation 99 we obtain 



u”! 1 = At 



l ' ] k 2 L h 2 



U n + U? 

U i+1,j + U i-\,j 
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Using the same method for Equation 98 we solve for v ^ to obtain 

.2 , „ a. 2 






a h 2 K w +v <- uV-ttHb' + 



(100) 
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( 101 ) 



To ensure the stability of Equations 100 and 101 the Von Neumann stability 
condition of 
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( 102 ) 



must be satisfied. 2 



D. APPLICATION OF THE STRESS-FREE BOUNDARY CONDITIONS TO 
THE SOLID 

The boundary conditions of the two-dimensional domain are broken 

( il'' 

down into two major categories, those whose normal vector is q , where 



0 



V" J 



x xx = Try = 0/ and those whose normal vector is |^ + |J where t = x X y = 0. 
These are in turn divided into two classes. For n = q they are 

m 

al. The boundary whose unit normal is q , that is facing in the positive x- 
direction, the surface El in Figure 5 and 



a2. The boundary whose unit normal is 
direction, the surface DH in Figure 5. 



(~1\ 

J 



, facing in the negative x 



Similarly for n = 



(0 ' 
±1 



( 0 \ 

bl. Boundary whose unit normal is ^ I , the surfaces AB, GH and IJ in 
Figure 5 and 



(0 \ 

b2. Boundary whose unit normal is , the surfaces CD, EF and KL in 
Figure 5. 

The application of the stress-free boundary conditions for cases al and bl is 
discussed below. 



2 For a brief treatment of the Von Neumann stability criteria see Appendix 
B. 
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1. Application of the Stress-free Boundary Condition for the Case of 




The boundary under investigation is identified in Figure 8 as XY. 
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Figure 8. Boundary with Normal 
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The governing equations as derived in Chapter II Section A 

du dv . 



(103) 
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(104) 
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(105) 



(106) 



Concentrating on displacements in the lateral (x) direction, a typical boundary 
point as shown in Figure 8 must satisfy the boundary conditions. Equations 
103 and 104 and the governing equations of motion. Equations 105 and 106. 

If we apply Equations 105 and 106 at the node (N,j) in Figure 8 we will 

require values for u n N+hj/ v n N+ i /y+1 ^N + l,H^N + l,;' w N + l,/ + l and 4+1,, -1' 
which lie outside of the domain and are called pseudonodes. To eliminate 
this reliance the technique as developed by Ilan and Lowenthal (Ilan, 1976, pp. 
431-453) is followed, and is presented here. 

The lateral displacement («) at the node in Figure 8 can be 

expanded in a Taylor series as 
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+ order / 
N ,j,n terms 



(107) 



where denotes the value of u at the node (N-l,j) and at time level n. 

Bu B 2 u 

Alternate expressions for ^ and are given by equations 104 and 105, we 
have from Equation 104 
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( 108 ) 
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From Equation 105, 
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( 110 ) 



Thus the lower order terms of Equation 107 can be written as 
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Using centered differences for and the following difference 



formula for the cross derivative term 

1 
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dxdy 2h 2 ^ ^0 + l N,j 1 N 1,;+1 
the finite difference approximation for Equation 111 is 
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Cancelling common factors Equation 113 reduces to 
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Solving for w N . explicitly we obtain 
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(115) 



The truncation error for Equation 115 is 0(/i 3 ) (Ilan, 1976, pp. 431-453). Using 

the same procedure for the transverse displacement an explicit expression 
n+1 . 

• is given as 
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2. Application of the Stress Free Boundary Condition for the Case of 




In the case of b2, the governing equations 



_ du dv 
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(119) 



( 120 ) 



and a typical boundary point is depicted in Figure 9. Expanding in the vertical 
direction at the node (i, M+l) and at time level n, and ignoring higher order 
terms 
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from Equation 119, Equation 120 now becomes 
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Figure 9. Boundary with Normal 
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We use central difference formulas for and for the cross 

d 2 V 



derivative term 



dxdy 
d 2 v 



the following finite difference formula is used 
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The analogous finite difference equation is now 
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Cancelling common factors yields 
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and solving for u ,• M explicitly, we obtain 
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Similarly for y ,• M we have 
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E. FINITE DIFFERENCE APPROXIMATIONS FOR THE CORNER NODES 

The I beam shaped domain has four 270° corners which are identified in 
Figure 10 as 1, 2, 3 and 4. The treatment of these corners falls into two 
categories, a) corners 1,3, and b) corners 2,4. Within each category the finite 
difference formula applied is identical, the difference between them comes in 
the cross derivative terms of the elastic wave equation. Each category will be 
discussed below. It is important to note that only the governing equations of 
motion are applied. The stress free boundary conditions are omitted due to 
the complexity that arises in trying to apply stress conditions at the corner 
node. We assume as in Fuyuki and Matsumoto that the consequences of 
neglecting these boundary conditions is minimal. (Fuyuki, 1980, pp. 2051- 
2069) 

Category a (Corners 1,3): An arbitrary numerical mesh with Ax = Ay = h 
about corner 1 is presented in Figure 11. The governing equations of motion 
are 
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Figure 11. Comer Node — Category a 



In Equation 124 which describes lateral motion, the normal central difference 
, d 2 u d 2 u d 2 u 

formulas are used for the / terms and for the cross derivative 



term 
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the difference formula of Fuyuki and Matsumoto (Fuyuki, 1980, 



pp. 2051-2069) is applied at the node (i,j) which is 
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Where D + is the forward difference formula 
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and D is the backward difference formula 
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The resultant difference formula for the cross derivative term in v at node i,j 
is 
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which is 0 (/j 2 ). The finite difference equivalent of Equation 124 is 
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Solving explicitly for from Equation 130 we obtain 
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Similarly for Equation 125 an explicit expression for v i j is 



44 



M + l . 

'■> kih 2 



5V*v)' 



2 - 



2AV 






1 1 






\ k l + k h) 






+ - 



At* 



2/d 



yti k T j 



L.n i J , ,,n o,,)i i „,n-l 

( u i+-[,j + u i-\ l j + u i / j+\ + u i,j-1 u i+l,j+1 2u i,j] v i,j ■ 



(138) 



For Category b), the governing equations of motion are the same. The finite 
difference approximation of all terms with the exception of the cross 
derivative are done in the same manner. For the cross derivative term in the 
case of Equation 124, the difference formula applied at the node (i,j) is 



dxdy 2 - 



D 



[D y + D*D y ] 



d 2 V 



Thus at that n ode in Figure 12 is represented by 
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The finite difference approximation to Equation 124 is now 
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Figure 12. Corner Node — Cateogory b 



n + 1 



Solving explicitly for u 2 - j we obtain 
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in the same manner an explicit expression for v 2 - j is 
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F. BOUNDARY CONDITIONS AT THE FLUID/SOLID INTERFACE 

The boundary conditions at the fluid/solid interface are modified by the 
introduction of a normal plane wave incident on the surface of the solid. As 
a result the normal component of stress is no longer zero. We must allow for 
the effect of the scattered pressure at the interface, which is a result of the 
compliance of the I-beam structure and reflections of displacement waves 
from internal boundaries. This is done by use of the compatibility condition 
as derived earlier. Our necessary equations are 
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(148) 



dp d 2 V 

dy dt 2 



2 ^P_^P + ^P 

f dt 2 dx 2 dy 2 



(149) 



where e = pf/ p s and p s is the scattered pressure along the boundary y = Kh (at 
the interface, see Figure 13), and at time level n. 

Explicit expression for and v are derived in the same manner as 
before and they are 
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Figure 13. Fluid-Solid Interface 



Looking at the compatibility condition. Equation 142 we use central difference 
approximations to obtain 



n n - n n 
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2 ft 
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(152) 



assuming that it is applied at the boundary y = Kh in Figure 13. 
Solving for we obtain 

„ 2ft 






Kk - lv i,K +V i,K J + Pi, 



K+l 



(153) 



all quantities on the right hand side are known. We now have a method of 
calculating which is required when applying Equation 149 at the 

boundary y = Kh. 

Solving for the p” +1 term of Equation 149 we have 
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The components of the right-hand side with the exception of lie in the 

interior or on the boundary of the fluid domain and can be evaluated, 
however p”^_ 1 is essentially a pseudonode for the fluid and is determined 

from Equation 153 above. By this method we have now generated the 
scattered pressure waves caused by the vibration of the solid at the fluid /solid 
surface. 



G. PROGRAMMING CONSIDERATIONS 

The program for this thesis was written entirely in Matlab 4.0 Beta 

version for two reasons, 

• Ease of programming 

• The ability to generate quality graphics. 

Although originally intended as a linear algebra toolkit the above features 
have caused it to be used more and more as a high level programming 
language. In our case Matlab was convenient since the fluid and solid 
domains are square matrices, which are easily manipulated in Matlab. 
Updating values is done somewhat differently than in FORTRAN or C and is 
discussed below. 

Equation 66 of this section updates the interior points of the fluid domain 
and is given by 




an equivalent FORTRAN statement might look like 
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DO 10 I = 2 , K 
DO 20 J = 2, K 



P ( I , J, N+l ) = 2.0* (1 .0-2 .0* (RHO**2) ) *P (I, J,N) 

& + (RHO**2) *(P(I-1,J,N) +P (1+1, J,N) +P ( I , J-l , N) +P ( I , J+l , N) ) 

& -P(I,J,N-1)) 

20 CONTINUE 
10 CONTINUE 

Here each value is updated individually by using a double DO loop. In 
Matlab this is done by "shifting" a grid the size of the interior around the 
appropriate matrix and weighting terms. The equivalent code in MATLAB 
would be 

PNEW (2:K,2:K) = 2*(l-2*RHO A 2)*PCURR(2:K,2:K) 

+(RHO A 2) It (PCURR(l:K-l,2:K)+PCURR(3:K+l,2:K)+PCURR(2:K,l:K-l) 

+PCURR(2:K,3:K+l))-POLD(2:K,2:K). 

where PNEW contains the new values, PCURR the current values and so on. 
All updating is done in this manner eliminating the requirement for 
multiple do loops. 
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IV. NUMERICAL RESULTS 



In an effort to verify our code we checked the behavior of the fluid and 
employed energy conservation methods to check for the consistency of the 
coupled domain. For the fluid waveguide we want to ensure that the 
propagating modes behave as expected, that is, they should not reflect from 
the artificial boundary. This was done by placing a driving force of the form 



Section A) at the boundary y = 0, which excited the fundamental, first and 
second modes (n = 0, 1, 2). The coefficients A n had value 1 for all n. The 
amplitude of each mode was measured at the boundary. As can be seen from 
Figures 14 a, b, and c the fundamental and first mode approach the value 1 
with the second mode showing the same behavior but at a much slower rate. 

To check that the coupling of the two domains was working correctly we 
eliminated the cavities of the I-beam which reduced our problem to one 
which could be solved analytically. For a normally incident plane wave only 
the fundamental mode was excited. Since the incident wave displaces the 
solid only in the y direction v has no x dependence and u is identically zero. 
The values of Aq and v are given by (Scandrett, 1992, interview). 



p= Aj^*rnX-t) 



(155) 




A) = C 1 




and 
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40 60 80 100 120 140 

Fig 14c time - seconds(scaled) 



For kf = 1 and kL = 0.2542, c\ and c 2 were given by 0.0625 + 0.03/ and 
0.0585 - 0.0374/ respectively. In this case the magnitude of ||Ao|| of 0.1211 
compared favorably with the numerical solution of 0.1255. A major 
discrepancy lay in the numerical and analytical values of v (transverse 
displacement). The average absolute value for the numerical solution was 
13.12 while the average absolute value for the analytical solution was 0.0965. 
They exhibited similar behavior but the numerical solution was translated by 
a constant term. We believe this to be a result of there being no displacement 
term to compensate for the effect of imposing a periodic pressure 
instantaneously in time which has caused our solid domain to drift or 
displace, violating one of our initial assumptions (see Chapter I, Section A). 
To nullify this effect we take the time derivative of our steady state solution 
for v and then compare with our analytic value as can be seen in Figure 15. 
Again the numerical and analytical solutions give close agreement. A second 
check was to apply energy conservation methods to our steady state solution, 
from which it can be shown (Scandrett, 1992) that the propagating modes 
must satisfy 



M q 1 

z/ywr =- 2 A) R e(A))=-- im \p 



n=-M 



-1 



y = 0 



V 



y = Qi 



dx. 



(156) 



2 

The quantities ^ n -_ M Pni\T anc * -2/?oRe(A 0 ) for the various values of kf 

are listed in Table 1. Considerable discrepancies exist throughout, which may 
indicate either an error in the code or the inability of our simulation to model 
the high frequency components which exist at the interface. Another factor 
which may contribute to the failure of the integral is the translation of v 
which was discussed previously. With this in mind we must possibly 
consider the following results as being inaccurate. 
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DOMAIN (DXS = 1/40) 



The first aspect we look at is the amplitude of the propagating modes of 
the scattered pressure. For kf = 7 , (this we call our original system) the time 
history of the amplitudes for the fundamental, first and second modes are 
shown in Figures 16, 17, and 18 respectively. These indicate roughly the rate 
of convergence of the numerical method. 



TABLE 1. ENERGY CONSIDERATIONS 



*/ 


( 1 \ 

1 1 

— Im f p ~v ~dx 

2 J/y = o y = o 
V-i / 


x1mA.ikii 2 


-2A,Re(A,) 


0.5065 


-0.1735 


0.022 


0.0003 


1 


0.475 


0.0279 


-0.2138 


1.2516 


0.12965 


0.0006 


0.055 


2.026 


40.5886 


2.135 


3.782 


3.3446 


-0.6093 


0.5380 


-1.3844 


4 


2.6288 


0.2179 


-0.6506 


4.5585 


7.6195 


0.2967 


0.5685 


6.5572 


-5.72 


0.0835 


-0.4053 


7 


4.0212 


0.0644 


0.0667 


7.5 


0.18715 


0.0846 


0.6723 


8 


4.13916 


0.5225 


0.2738 


9 


4.8779 


0.0742 


-1.5974 



To gain insight into the characteristics of the solid, we treated the flange at 
the fluid-solid interface as a thin plate and tried different values of kf 
corresponding to the resonant frequencies of a thin plate with prescribed 
boundary conditions. The first case was a plate with the fixed (referred to as 
FF) boundary conditions given by 

v{-\) = i;(l) = i?"(-l) = v"(l) = 0. 
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Amplitude of the fundamental mode 
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Amplitude of the first mode 
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Amplitude of the second mode 
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Fig 18 time - seconds(scaled) 



The resonant frequencies are 



cjhn 2 n 2 

| — - 

n 4a^6(l- u ) ' 



(157) 



(Lalanne, 1982, p. 103) where v is Poisson's ratio, h is the thickness of the 
plate, cj is the transverse velocity of the solid and a is a scaling constant for 
distance. 

For the second case, clamped (referred to as CC) boundary conditions were 
prescribed, which are given by 

i>(-l) = v{\) = i>'(-l) = v'{\) = 0. 



The resonant frequencies of this system are given by 

0, __E7 

" 4^6(1 - v) 



(158) 



where X? = 22.37, xf = 61.67 and xf = 120.9 . (Lalanne, 1982, p. 103) 

The purpose of this experiment is to determine whether the I-beam 
exhibits behavior similar to a clamped or fixed plate, since it would be 
advantageous (numerically) if we could substitute a periodically placed 
boundary condition for the remaining portion (lower flange and center spar) 
of the I-beam rather than having to calculate finite difference approximations 
for the entire I-beam. 

We do this by plotting the amplitudes versus the corresponding values of 
kf for the fundamental, first, and second modes. The presence of any peaks 
would indicate a resonant type behavior. If any of these peaks corresponds to 
a particular value of kf for a clamped or fixed plate we say that for that value 
of kf (and hence &>), the I-beam behaves in a manner similar to a plate with 
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clamped or fixed boundary conditions. This gives us a possible range for kf 
values on which to concentrate when looking for resonant frequencies. 

The values of kf are given in Table 2. The quantities in parentheses are 
the modes that propagated for a particular value of kf which is determined by 
the maximum integer value n can take such that (3 n is real (Ref. Equation 71). 



TABLE 2. k f VALUES 





SYSTEM 


CATEGORY 


ORIGINAL 


c-c 


F-F 


1 


1(0) 


1.2516(0) 


0.5265 (0) 


2 


4 (0, ±1) 


3.3446 (0, ±1) 


2.026 (0) 


3 


7 (0, ±1, ±2) 


6.5572 (0, ±1, ±2) 


4.558 (0, ±1) 



We divide our values of kf into three categories (for identification 
purposes). In the first category we include the values of kf corresponding to 
the first resonant frequency of the clamped and fixed plates and the values of 
kf for our original system which allows only the fundamental mode to 
propagate. In the second category are values of kf that correspond to the 
second resonant frequency of the clamped and fixed plates as well as the value 
of kf which allow only the fundamental and first modes to propagate. The 
third category corresponds to the third resonance of the clamped and fixed 
plates and the first three modes in our original system. Also included in this 
final category are the values of kf = 7.5, 8 and 9 which will allow us to study 
the behavior of the second mode at higher frequencies in greater detail. 
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In Figure 19 we show a plot of kf values versus the amplitude of the 
fundamental mode. Computed values are shown in black. A cubic spline of 
the points involved is shown in red. We must note however that this spline 
is only a possible representation of the behavior of the amplitude for different 
values of kf, since time constraints did not allow us to conduct a more 
comprehensive set of simulations from which we could obtain an accurate 
picture of the behavior. 

It can be seen that for kf = 0.5065 and 2.026 (the points labeled FF1 and FF2, 
the I-beam is exhibiting a resonant type behavior. These values correspond to 
the first and second resonant frequencies of a fixed plate. It can also be seen 
that for values of kf greater than 3.5 there is little variation in the amplitude 
of the fundamental model since we are now past the first cutoff value (after 
which three modes propagate) of kf = k. This leads us to believe that at these 
higher frequencies most of the activity takes place in the higher modes. 
However the total energy calculation shifts from fundamental to first and 
back to fundamental as can be seen in Table 3. 

In Figure 20 the only resonant behavior is exhibited for the value 
kf= 3.446 (the point labelled CC2). This frequency is just past the first cutoff 
value (n) and the amplitude of the first mode is double that of the 
fundamental (compared with the corresponding value in Figure 19) which 
confirms that at higher frequencies energy is being propagated in the higher 
modes. 

In Figure 21 we plot the amplitude of the second mode against kf. The 
only resonant type behavior which exists here is for the kf = 8, but since the 
amplitudes involved are so small it would be difficult to draw an accurate 
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kf vs Amplitude of fundamental mode 
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Fig 19 kf values 
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Fig 20 kf values 



kf vs Amplitude of second mode 
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Fig 21 kf values 



conclusion about the existence of a resonant frequency. Another aspect we 
investigate is the proportion of the total energy carried by each propagating 
mode which is given by the formula 




p=-M 



(159) 



where E n is the energy of the n th mode and A n the amplitude of the n th mode 
(Kinsler, 1982, p. 110). These quantities are summarized in Table 3. 



TABLE 3. ENERGY CONSIDERATIONS 





£„(%) 


kf 


£ o 


Ei 


E 2 


£-i 


£-2 


0.5056 


100 


N/A 


N/A 


N/A 


N/A 


1 


100 


N/A 


N/A 


N/A 


N/A 


1.2516 


100 


N/A 


N/A 


N/A 


N/A 


2.026 


100 


N/A 


N/A 


N/A 


N/A 


3.3446 


13.04 


43.47 


N/A 


43.47 


N/A 


4 


13.36 


40.82 


N/A 


40.82 


N/A 


4.5535 


7.47 


46.26 


N/A 


46.26 


N/A 


6.5572 


82.43 


2.29 


6.5 


2.29 


6.5 


7 


58.23 


8.08 


12.81 


8.08 


12.81 
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We can see that the proportion of energy carried by the fundamental mode 
varies from 100% to 12.04% as kf goes from 2.026 to 3.3446. This can be 
accounted for as a redistribution of energy that takes place as the frequency of 
the waves passing through the fluid wave guide pass the first cutoff value of 
kf = k. Similar arguments hold for higher values of kf. 

Finally we take a brief look at the shear strain field generated in the solid. 
We do this for two reasons, one to check the behavior at the corner nodes and 
two, to find out where the maximum shear strain occurs. 

When treating the corner points we neglected applying the stress free 
boundary conditions there assuming this would have no adverse effect on 
the behavior of the solid. Were this assumption invalid we would expect to 
observe singularities or other irregular behavior. Provided in Figures 22, 23 
and 24 are snapshots of the solid at different time levels which display no 
unusual behavior at the corner nodes, leading us to accept the assumption as 
valid. 

From an engineering standpoint we are interested in where the 
maximum shear strain occurs for possible failure analysis. As can be seen 
from Figure 25 the maximum occurs along the transverse borders of the 
I-beam cavity. We must note that the amplitudes and frequencies of the 
acoustic pressure waves involved in our study are not comparable to those 
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which would cause permanent deformation or failure, but which could result 
in fatigue cracking if subjected to prolonged periods of stress. 
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Shear Strain at time t = 3s(scaled) 
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Fig 23 Shear Strain at time t = 6s(scaled) 
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Fig 24 Shear Strain at time t = 12s(scaled) 
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hear Strain at steady state 



V. CONCLUSIONS 



Several problems arose in implementing the numerical scheme. It was 
found that the non-local radiation boundary condition did not completely 
annihilate the initial wavefront and a small amount of reflection took place, 
but over time the effect of this reflection became negligible. Evanescent 
modes did not decay sufficiently and were reflected at the boundary adding to 
the overall noise of the problem. 

When calculating the amplitudes of the propagating modes at steady state 
we applied the formula 



which should yield consistent results for any y values in the fluid domain. In 
the neighborhood of the fluid/solid interface this was not the case as can be 
seen from Figure 26. We believe this is due to the presence of high frequency 
components in this region and having too coarse a mesh to effectively 
evaluate the integral there. 

With a stepsize h of 1/40 our truncation error was on the order of 10 -3 . 
To increase the accuracy we can perform Romberg extrapolation or decrease 
the stepsize, thereby increasing the dimensions of the matrices involved. The 
latter was not an option due to the construction of the code and machine 
limitations, in that a simulation which involved 10 to 15 thousand time steps 
took from 10 to 12 hours to perform. Increasing the size of the matrices 




l 



(161) 



-1 
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involved, thereby increasing the required time, would not make it a suitable 
code for experimentation and timely results. 



Amplitude of fundamental mode 




Figure 26 . Amplitude of Propagating Modes in the Y-domain 

Aspects of the problem which deserve further study would be the use of 
obliquely incident waves vice normal as was used here. Variation of the 
cavity size and its effect on the amplitude of the propagating modes should 
also be considered. As a simplifying assumption the cavities were considered 
to be void. It would be realistic to expect that they may contain fluid 
(seawater) or an acoustic dampening material. Modifying the model to 
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account for such conditions warrants further study. A final aspect that could 
be considered is a resonance analysis. 

Throughout all of our finite difference approximations we were able to 
maintain second order accuracy in space and time. We did not have to resort 
to the use of pseudonodes outside of our fluid-solid domain. Finally the 
possibility of being able to establish an acoustic signature for a double-hulled 
structure could open up a new avenue of submarine detection for which this 
thesis could be considered a starting point. 
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APPENDIX A. VON NEUMANN STABILITY ANALYSIS FOR THE 2-D 



SCALAR WAVE EQUATION 



The general form of the scalar wave equation being used is 

1 B 2 u d 2 u d 2 u 



c) dt 2 Bx A By' 



(A-l) 



whose equivalent difference equation, given Ax = Ay = h is 

1 ■("" J 1 - Ki + ""7') - + <v)+£K-* i - Hr 

(A-2) 



cjAt 2 



The error function takes the form 

u = E$ y'(0p + w)\ 



(A-3) 



This is substituted into Equation A-2 and common terms are cancelled to give 
^(S-2H-') = ±(e* h -2 + e-W)^(e‘*-2 + e->*). (A-4) 



Using the following identity. 



cos OCX = 



e iax +e -iax 



(A-5) 



letting ph = and defining p = cAt/h, Equation A-4 reduces to 

£ -2+ £,~ x = p 2 (Acosfih-A). 

Multiplying Equation A-6 through by £ and collecting like terms gives 



(A-6) 
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£ 2 -2^1-2p 2 (cosj3p-l)j + l = 0 



(A- 7) 



or equivalently 



£ 2 -2<^l-4p 2 sin 2 ^^ 



+ 1 = 0 . 



The roots of this quadratic are 



^1^2=1- 4p 2 sin ^y] ± ~ 4 P 2 sin2 (y 



v \2 



For stability we require 



which forces 






i,2. 2 (PhX) 

l-4p sin | — 






- 4 < 0. 



Solving for p 2 in Equation A-10 gives 



P 2 <— 1 
2 sin 2 



ph 

U. 



which reduces to our final stability criterion of 

p. 1 



J2sin 2 f— 



or 






V2 



(A-8) 



4 . (A-9) 



(A-10) 



(A-ll) 



(A-12) 
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APPENDIX B. VON NEUMANN STABILITY ANALYSIS FOR THE ELASTIC 

WAVE EQUATION 



In this section we give a brief outline of the stability analysis required for 
the elastic wave equation, whose general form is 



2 d 2 U 2 <^"U I 2 2 \ 3 2 W 

CL—J + Ct—J + (Cl-Ct) 
dx dy v ' 



dxdy dt 2 



(B-l) 



9 ? 9 9 

2 d V 2 O V (2 2 \ oU d V 

c r — + C L — + f 



dx z 



3/ 



(cl-ct) 



dxdy dt 2 



The finite difference approximation for Equation B-l is 



(B-2) 



%( u M.r 2u "i +u i-u) + tf( u ”H - 2 ""/ + “" h ) 

/ 2 _ 2 \ 

~ v m.h - y+i + = ^2 (""/ 1 - 2 “"; + <B ‘ 3) 



We use the following error functions 

u = US r e i (P p+ w) h and v=V^ r e i ^ p+ ^h (B-4) 

which are substituted into Equation B-3, common terms are cancelled and 
complex exponentials are gathered into trigonometric quantities to give 

-4 r 2 (cl sin 2 y + sin 2 &-\l - r 2 (cl - c£)(sin £/zsin yh)V = ($ - 2 + S' 1 ) • U.(B-5) 



Following the same procedure for Equation B-2 we obtain 
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\ 






fin 


= A 


fin 








/ 







Equations B-5 and B-6 can be written in matrix form which is 

-4 r 2 ^ c 2 sin 2 ^- + Cj sin 2 j -r 2 {cl - c 2 jsin (3h sin yh 
-r 2 {cl -cj'jsinph sin yh -4r 2 sin 2 ^ + cl sin 2 -y- 



where A = t, - 2 + The eigenvalues of matrix in Equation B-7 are 

r 2L2_ c 2\ 

2 r 2 {cl + cj S j(cosph + cosyii-2 )± — — — ^( cos(/?/i - yh) + cos(ph+ yh)-2 ) 



they take on a maximum value when ph = yh = n, and we obtain 

At = ^2 = -4r 2 ^c 2 + cj ). 

We are now left with the identity 

^-2 + r 1 =-4r 2 (c£+c?) 



(B-7). 



1 

2 ^2 



(B-8). 



(B-9) 



(B-10) 



or 



£ 2 + ^4r 2 |c 2 +c 2 j- 2j^ + 1 = 0. 



For stability we require 



(B-ll) 






which forces 



or 



^4r 2 (c^ + cj')-2j -4 < 0 



(B-12) 



At < 



h 



(B-13) 
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APPENDIX C FINITE DIFFERENCE FORMULAE FOR THE EQUATION OF 
MOTION AND BOUNDARY CONDITIONS 



This appendix contains the finite difference approximations of all the 
equations required for the simulation. 

Solid 

1) Elastic Wave Equation 
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Stress-free Boundary Conditions 
Surfaces with n = fj 1 
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Surfaces with n = L 
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Corners 2 and 4 
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Fluid 

1) Wave Equation 




Radiation Boundary Condition 
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APPENDIX D. COMPUTER CODE 



% function basel - base for the ' I - beam ' shaped domain 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% constructs a grid in the shape of an I beam 

function [b,rows, pcoll,pcolr, fill , coll , coir] = basel(n) 

% variables 

% rows - identifies the rows to be zeroed to form the I beam 

% coll - identifies the columns to be zeroed to the left of the center spar 

% coir - identifies the columns to be zeroed to the right of the center spar 

% pcolr - identifies the columns to be zeroed to the right of the center spar 
% including the boundary values. 

% pcoll - identifies the columns to be zeroed to the left of the center spar 
% including the boundary values. 

% bb - building block of the correct size 

% cc - dimension of the blocks to be zeroed out on either side of the center 
% spar 

% s - variable to adjust the size of the zero blocks and keep it symmetric 
% fill - block of zeros of appropiate size to fill in the spaces on 

% either side of center spar i.e. zeroing out, at every time step. 

% m - no. of divisions in the half length of the domain. 

% cr - variable used to pick the columns to the right of the spar 



% build a grid of the correct size, 
bb = ones (2*n+l) ; 



% the basic variables required to build the I - beam shaped grid 
m = n+1; 
s = log2 (n) ; 
cc = (s-1) * (log2 (n) ) ; 

% determine the row numbers to be zeroed out. 
rows= m- (cc-1) :m+ (cc-1) ; 



% columns to the left of the center spar, pcoll includes the boundary points 
% coll does not. 
coll = l:cc+l; 
pcoll = l:cc+2; 



% columns to the right of the center spar, pcolr includes the boundary points 
% coir does not. 

cr = 2 *n+l - (cc) ; 
coir = cr:2*n+l; 
pcolr = cr-l:2*n+l; 

% construct the block of zeros used to fill in the spaces on either side 
% of the center spar at every time step, 
e = size(rows); e = e(2); 
d = size(coll); d = d(2); 
fill = zeros (e , cc+1 ) ; 
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bb(rows,coll) = fill; 
bb (rows , coir) « fill; 

b = bb; 
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% function bndfnx - boundary facing in the negative x direction 
% i.e unit normal = (-1,0) 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% bndfnx applies the traction free boundary and the elastic 
% wave equation along the boundary of the I-beam with unit normal 
% (-1,0) the technique used is that as developed by Ilan and Lowenthal 

% and is not discussed here. 

% the columns containg the nodes on the surface in question and 
% the required neighbours (those one column in from the surface) 

% are picked off from the current and old values of u and v 
% the shifted and weighted with constants from the vectors 
% cunx (Constants for U-values for the Negative X boundary) 

% and cvnx (Constants for V-values for the Negative X boundary) 

% and inserted in the correct position of the updated u and v. 



function [un ,vn] = bndf nx (uc , vc , uo, vo , un , vn, rows , pcoll , c, d) ; 



% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcoll carries out the same function as rows for 
% the columns, and the contains the column location of the 
% of the boundary facing the negative x direction 

% c = cunx; 

% d = cvnx; 



[ i3 , j 3 ] = size(rows); 

mrows = [rows(i3)-l rows rows ( j 3 ) +1 ] ; 

[i4,j4] = size(mrows); 

[i5,j5] = size(pcoll); 

% cul cu2 cvl cv2 co cov contain the necessary u and v values 
% for our calculations 

cu2 = uc (mrows , pcoll ( j 5) +1 ) ; 
cul = uc (mrows , pcoll ( j 5) ) ; 

cv2 = vc (mrows, pcoll ( j5) +1) ; 
cvl = vc (mrows , pcoll ( j 5) ) ; 

co = uo (mrows , pcoll ( j5) ) ; 
cov = vo (mrows , pcoll ( j 5)); 



% the updated values are calculated 

uc 1 = c ( 1) *cul (2 : j4— 1) - co (2 : j4-l) + c ( 2 ) *cu2 ( 2 : j 4-1) . . . 

+c ( 3 ) * (cul ( 3 : j 4 ) + cul (1: j4-2) ) . . . 

+c ( 4 ) * ( cv2 ( 1 : j 4 -2 ) - cv2 ( 3 : j 4 ) ) . . . 



88 



+c(5)*(cvl(3: j4) “ cvl ( 1 : j4-2) ) ; 



vcl = d ( 1 ) *cvl ( 2 : j 4-1 ) - cov ( 2 : j4-l) + d (2 ) *cv2 (2 : j4-l) . . . 
+d ( 3 ) * ( cvl ( 3 : j 4 ) + cvl (1 : j4-2) ) . . . 

+d(4) * (cu2 (l: j4-2) - cu2(3: j4) ) . . . 

+d (5) * (cul ( 3 : j 4 ) - cul (1: j4-2) ) ; 



% and put in their proper place in un and vn 

pb= size (mrows) ; 

pbc = mrows (2 : pb (2 ) -1 ) ; 



un (pbc , pcoll ( j5) ) = ucl; 
vn (pbc, pcoll ( j5) ) = vcl; 
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% function bndfny - boundary facing in the negative y direction 
% i.e unit normal = (0,-1) 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% bndfny applies the traction free boundary and the elastic 
% wave equation along the boundary of the I-beam with unit normal 
% (0,-1) the technique used is that as developed by Ilan and Lowenthal 

% and is not discussed here. 

% the rows containg the nodes on the surface in question and 
% the required neighbours (those one row in from the surface) 

% are picked off from the current and old values of u and v 
% the shifted and weighted with constants from the vectors 
% cuny (Constants for U-values for the Negative Y boundary) 

% and cvny (Constants for V-values for the Negative Y boundary) 

% and inserted in the correct position of the updated u and v. 



function [un ,vn] = bndfny (uc,vc,uo,vo,un, vn, rows, pcoll,pcolr, c,d) ; 

% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcoll and pclor are both required as there are two regions, one on 
% either side of the center spar of the I -beam which require 
% our attention and they contain the location of the nodes in question 



% c = cuny; 
% d = cvny; 



[il,jl] = size(pcoll); 

[ i2 , j 2 3 = size (rows); 

[sr sc] = size(uc); 

% rl** ru* rv* and ro* pick off the rows on either side of the 
% center spar of the necessary u and v values. 

rlul = uc (rows (j2) +1 , pcoll) ; 
rlu2 = uc (rows ( j2) +2 , pcoll) ; 
ru2 = uc (rows (j2 ) +2 , pcolr) ; 
rul = uc (rows (j2 ) +1 , pcolr) ; 

rlvl = vc (rows ( j2 ) +1 , pcoll) ; 
rlv2 = vc (rows ( j2 ) +2 , pcoll) ; 
rv2 = vc(rows (j2)+2 , pcolr) ; 
rvl = vc (rows (j2) +1, pcolr) ; 

rol = uo(rows(j2)+l, pcoll) ; 
ro = uo(rows (j2) +1, pcolr) ; 
rolv = vo(rows(j2)+l, pcoll) ; 
rov = vo (rows (j2) +1, pcolr) ; 
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% the updated values are calculated 



ul = c(l) *rlul(2: jl-1) - rol(2:jl-l) + c (2) *rlu2 (2 : j 1-1) . . . 
+c (3 ) * (rlul (3 : j 1) + rlul ( 1: j 1-2) ).. . 

+c (4 ) * (rlv2 ( 1 : j 1-2 ) - r lv2 ( 3 : j 1) ) . . . 

+c (5) * (rlvl (3 : j 1) - rlvl (1: jl-2) ) ; 

ur = c(l) *rul (2: jl-1) - ro(2:jl-l) + c (2) *ru2 (2 : jl-1) . . . 

+c ( 3 ) * ( rul ( 3 : j 1 ) + rul (1: jl-2) ) . . . 

+c(4) * (rv2 (l: jl-2) - rv2 (3 : j 1) ) . . . 

+c (5) * (rvl (3 : jl) - rvl ( 1 : j 1-2 ) ) ; 



vl = d ( 1) *rlvl(2: jl-1) - rolv(2:jl-l) + d ( 2 ) *rlv2 ( 2 : jl-1) . . . 
+d(3)*(rlvl(3:jl) + rlvl(l: jl-2) ) . . . 

+d(4) * (rlu2(l: jl-2) - rlu2 (3: jl) ) . . . 

+d (5) *(rlul(3:jl) - rlul(l: jl-2) ) ; 

vr = d(l)*rvl(2: jl-1) - rov(2:jl-l) + d (2) *rv2 (2 : jl-1) . . . 

+d ( 3 ) * (rvl ( 3 : j 1 ) + rvl ( 1 : j 1-2 ) ) . . . 

+d(4) * (ru2 (1: jl-2) - ru2 (3 : j 1) ) . . . 

+d(5) * (rul (3 : jl) - rul ( 1 : jl-2)); 



% and put in their proper place in un and vn 

pb= size(pcoll) ; 

pbl = pcoll ( 2 : pb ( 2 ) -1) ; 

pbr = pcolr (2 :pb(2) -1) ; 



un (rows ( j 2 ) +1 , pbl ) = ul; 
un (rows ( j2) +1 ,pbr) = ur; 



vn (rows ( j2 ) +1 , pbl) = vl; 
vn (rows ( j 2 ) +1 , pbr) = vr; 

% the same procedure is reppeated for the 'bottom' of the I-beam 

un ( 1 , 2 : sc-1 ) = c ( 1) *uc ( 1 , 2 : sc-1) - uo ( 1 , 2 : sc-1) . . . 

+ c ( 2 ) *uc (2,2: sc-1) . . . 

+c (3) * (uc ( 1, 3 : sc) + uc ( 1 , 1 : sc- 2 ) ) . . . 

+c ( 4 ) * ( vc ( 2 , 1 : sc-2 ) - vc ( 2 , 3 : sc) ) • . . 

+c (5) * (vc ( 1 , 3 : sc) - vc (1, 1: sc-2) ) ; 



vn ( 1 , 2 : sc-1 ) = d (1) *vc (1, 2 : sc-1) - vo ( 1 , 2 : sc-1) . . . 
+ d(2) *vc(2,2:sc-l) . . . 

+d (3) * (vc (1 , 3 : sc) + vc (1 , 1 : sc-2) ) . . . 

+d ( 4 ) * (uc ( 2 , 1 : sc-2 ) - uc (2 , 3 : sc) ) . . . 

+d (5) *(uc(l / 3:sc) - uc ( 1 , 1 : sc-2) ) ; 
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% function bndfpx - boundary facing in the positive x direction 
% i.e unit normal = (1/0) 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% bndfnx applies the traction free boundary and the elastic 
% wave equation along the boundary of the I-beam with unit normal 
% (0,1) the technique used is that as developed by Ilan and Lowenthal 

% and is not discussed here. 

% the columns containg the nodes on the surface in question and 
% the required neighbours (those one column in from the surface) 

% are picked off from the current and old values of u and v 
% the shifted and weighted with constants from the vectors 
% cupx (Constants for U-values for the Positive X boundary) 

% and cvpx (Constants for V-values for the Positive X boundary) 

% and inserted in the correct position of the updated u and v. 

function [un ,vn] = bndfpx(uc, vc,uo, vo,un, vn, rows,pcolr , c,d) 

% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcolr carries out the same function as rows for 
% the columns, and the contains the column location of the 
% of the boundary facing the positive x direction 



% c = cupx; 

% d = cvpx; 



[ i3 , j 3 ] = size(rows); 

mrows = [rows(i3)-l rows rows(j3)+l]; 
( i 4 , j 4 ] = size(mrows); 



% cul cu2 cvl cv2 co cov contain the necessary u and v values 
% for our calculations 



cu2 = uc (mrows , pcolr ( 1) -1 ) ; 
cul = uc (mrows , pcolr ( 1) ) ; 

cv2 = vc (mrows, pcolr ( 1) -1) ; 
cvl = vc (mrows, pcolr (1) ) ; 

co = uo(mrows, pcolr ( 1) ) ; 
cov = vo (mrows, pcolr (1) ) ; 

% the updated values are calculated 



ucr = c(l) *cul(2: j4-l) - co(2:j4-l) + c (2 ) *cu2 ( 2 : j4-l) . . . 
+c(3) * (cul (3 : j4 ) + cul (1: j4-2) ) . . . 

+c ( 4 ) * (cv2 ( 1 : j 4 -2 ) - cv2 ( 3 : j 4 ) ) . . . 

+c (5) * (cvl (3 : j4 ) - cvl (1: j4-2) ) ; 
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vcr = d ( 1) *cvl (2 : j 4 -1) - cov(2:j4-l) + d (2 ) *cv2 (2 : j4-l) . . . 
+d (3 ) * (cvl ( 3 : j 4 ) + cvl(l: j4-2) ) . . . 

+d (4 ) * (cu2 ( 1 : j4-2) - cu2 ( 3 : j4 ) ) . . . 

+d (5) * (cul (3 : j4 ) - cul ( 1 : j4-2 ) ) ; 



% and put in their proper place in un and vn 

pb= size (mrows) ; 

pbc = mrows (2 :pb(2) -1) ; 



un(pbc, pcolr ( 1) ) = ucr; 
vn (pbc, pcolr ( 1) ) = vcr; 
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% function bndfny - boundary facing in the positive y direction 
% i.e unit normal = (0,1) 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% bndfpy applies the traction free boundary and the elastic 
% wave equation along the boundary of the I-beam with unit normal 
% (0,1) the technique used is that as developed by Ilan and Lowenthal 

% and is not discussed here. 

% the rows containg the nodes on the surface in question and 
% the required neighbours (those one row in from the surface) 

% are picked off from the current and old values of u and v 
% the shifted and weighted with constants from the vectors 
% cupy (Constants for U-values for the Positive Y boundary) 

% and cvpy (Constants for V-values for the Positive Y boundary) 

% and inserted in the correct position of the updated u and v. 

function [un ,vn] = bndfpy (uc, vc,uo, vo,un, vn, rows, pcoll,pcolr, c,d) ; 

% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcoll and pclor are both required as there are two regions, one on 
% either side of the center spar of the I -beam which require 
% our attention and they contain the location of the nodes in question 



% c = cupy; 

% d = cvpy; 



[il,jl] = size(pcoll); 

[ i2 , j 2 ] = size(rows); 

% rl** ru* rv* and ro* pick off the rows on either side of the 
% center spar of the necessary u and v values. 



rlul = uc(rows (l) -l, pcoll) ; 
rlu2 = uc(rows ( 1) -2 , pcoll) ; 
ru2 = uc (rows ( 1) -2 , pcolr ) ; 
rul = ucjrows (1) -l , pcolr ) ; 



rlvl = vc(rows (1) -l,pcoll) ; 
r lv2 = vc (rows ( 1) -2 , pcoll) ; 
rv2 = vc (rows (1) -2 , pcolr) ; 
rvl = vc (rows ( 1) -1 , pcolr ) ; 

rol = uo(rows(l) -1, pcoll) ; 
ro = uo (rows ( 1) -l , pcolr) ; 
rolv = vo (rows ( 1) -1 , pcoll) ; 
rov = vo (rows ( 1) -1 , pcolr) ; 
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% the updated values are calculated 



ul = c ( 1) *rlul ( 2 : jl-1) - rol ( 2 : jl-1) + c(2) *rlu2 (2 : jl-1) . 
+c (3 ) * (rlul (3 : j 1) + rlul ( 1 : j 1-2 )).. . 

+c (4) * (rlv2 (l:jl — 2) - r lv2 (3 : j 1 ) ) . . . 

+c(5) * (rlvl (3 : jl) - rlvl (1 : j 1-2 ) ) ; 

ur = c(l)*rul(2: jl-1) - ro(2:jl-l) + c (2) *ru2 (2 : jl-1) . . . 
+c(3)*(rul(3: jl) + rul (1: jl-2) ).. . 

+c(4)*(rv2(l: jl-2) - rv2 (3 : jl) ) . . . 

+c ( 5 ) * ( rvl ( 3 : j 1 ) - rvl ( 1 : j 1-2 ) ) ; 



vl = d(l)*rlvl(2: jl-1) - rolv(2:jl-l) + d (2 ) *rlv2 ( 2 : j 1-1) 
+d(3)*(rlvl(3: jl) + rlvl ( 1 : jl-2 )).. . 

+d(4) *(rlu2(l: jl-2) - rlu2 (3 : jl) ) . . . 

+d (5)*(rlui(3:jl)- rlul (1: jl-2) ) ; 

vr = d (1) *rvl (2 : jl-1) - rov(2:jl-l) + d (2 ) *rv2 (2 : jl-1) . . . 
+d ( 3 ) * (rvl ( 3 : j 1) + rvl (1: jl-2) ).. . 

+d(4) *(ru2(l: jl-2) - ru2 (3 : jl) ) . . . 

+d (5) * (rul (3 : jl) - rul ( 1 : jl-2)); 



% and put in their proper place in un and vn 

pb= size (pcoll) ; 

pbl = pcoll (2 :pb(2) -1) ; 

pbr = pcolr (2 :pb (2) -1) ; 



un (rows ( 1) -1 , pbl) = ul; 
un (rows (1) -1, pbr) = ur; 

vn (rows ( 1) -1 , pbl ) = vl; 
vn (rows (1) -1 , pbr ) = vr; 
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% function craodck - checks the amplitude of the propagating mode 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% emodek is the driver program for the problem. All the vaue of the constants 
% are defined here and all quantities are scaled before they are fed into the 

% decks are cleared before calculation 
record ( ' erase ' ) 
clear 
clg 

axis ( 'auto' ) 

dx = input ('step size = '); 

a = input (' length scaling factor = '); 

omg = input (' time scaling factor (omega) = '); 

% constants 

% ct: tranverse velocity of the solid (steel) 

% cl: longitudinal velocity of the solid (steel) 

% dnss: density of the solid (steel) 

% cf : speed of sound in fluid (seawater) 

% dnsf : density of fluid (seawater) 

% epss: ratio ot fluid to solid densities 

ct = 3200; cl = 5900; cf = 1500; dnsf = 1026; dnss « 7700; 
epss = dnsf /dnss; 

% the determination of the scaled variables 
% dxs scaled distance 
% dts scaled time 
dxs = dx/a; 

kf = omg/cf;kt = (omg*a)/ct; kl = (omg*a)/cl; 

% dtsl and dts2 are the stability criteria for the solid and fluid 

% always choose the minimum 

dtsl = (kl*dxs) / (sqrt ( 1 + (ct*2/cl~2) ) ) ; 

dts2 = . 5* ( (kf *a*dxs) /sqrt (2 ) ) ; 

[dtsl dts2 ] 

dts = input (' desired time step = '); 
nn = input (' number of timesteps = ' ); 

% st - when to stop building the radiation boundary condition 
% for in the construction of the matrix A each loop thru the construction 
% cylce allows another mode to propagate 
st = stop(kf) ; 

i = sqrt (-1) ; 

% x a vector the length of the domain used in several places 
% i.e. when integrating etc. 

% mm - the mode being checked 0-fundamental etc. 

% nt - no of intervals in domain 
% tt weighting factor for the trapezoidal rule 
x = -1 : dxs : 1 ; 

m = size(x) ; m = m(2) ; ml = zeros (m) ; 

mm = 0 

nt = (1/dxs) ; 
tt = (m-1 ) / 2 ; 
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[un, vn , inn, zO ] = cw4 (ml , kl , kt , kf , dts , dxs ,x,nn,st, epss , tt , nt , mm) ; 



% un displacment in the x dirn 
% vn displacment in the y dirn 
% mn displacment of the fluid 

% zO vector containing the amplitude of the propagating mode 
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% function coupl - couples the two media at the fluid solid interface 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% coupl applies the traction boundary conditions at the fluid solid 
% interface where the shear component is zero but the normal component is 
% given by tau(yy) = -p (total) , which causes the I beam to deform 
% and causing the propagation of scattered pressure waves to propagate in 
% fluid. To prevent cavitation at the interface we apply the inviscid 
% form of the Navier-Stokes equation ( or Euler's equation ) at the 
% the boundary 



function [un ,vn,mnpd] = coupl (uc,vc,uo,vo,un, vn , c , d , k , dts , dxs , all , bll , epss ) ; 
% variables 

% un : updated values of u vn : updated values of v 
% uc : current values of u vc : current values of v 
% uo : old values of u vo : old values of v 



% c = cupy; 
% d = cvpy; 



i = sqrt (-1) ; 

[sr sc] = size (uc) ; 

% the forcing function 

time = exp(-i*k*dts) *ones ( 1 , sc) ; 




uc (sr, 1 : sc-2 ) ) . . . 

~ ) - vc (sr-1 , 3 : sc) ) < 
c(sr,l: sc-2 ) ) ; 



% the periodic boundary condition for u 
un(sr,l) = c ( 1) *uc (sr , 1) - uo(sr,l) . . . 

+ c ( 2) *uc (sr-1 , 1) . . . 

+ c ( 3 ) * ( uc ( sr , 2 ) + uc (sr , sc-1) ) . . . 

+ c ( 4 ) * ( vc ( sr-1 , sc-1 ) - vc (sr-l , 2) ) . . . 
+ c ( 5) * ( vc (sr , 2 ) - vc (sr , sc-1) ) ; 



un(sr,sc) = un(sr,l) ; 



% the normal component of stress plus the incident , ref lected 
% and scattered pressures 

vn (sr , 2 : sc-1) = d ( 1 ) * vc ( sr , 2 : sc-1 ) - vo (sr , 2 : sc-1) . . . 

+ d (2 ) *vc (sr-1 , 2 : sc-1) . . . 

+d (3) * (vc (sr , 3 : sc) + vc (sr , 1 : sc-2) ) . . . 

+d (4) * (uc (sr-1, 1: sc-2) - uc (sr-1 , 3 : sc) ) . . . 

+d ( 5) * (uc ( sr , 3 : sc) - uc (sr , 1 : sc-2) ) . . . 

+ ( (2*dts A 2 ) /dxs) * (2*epss*time (2 : sc-1) + epss*al 1 ( 2 : sc-1 ) ) ; 



% the periodic boundary condition for v 
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vn(sr,l) = d(l) *vc(sr, 1) - vo(sr,l) ... 

+ d ( 2 ) *vc (sr-1 , 1) . . . 

+d(3) * (vc(sr,2) + vc (sr , sc-1) ) . . . 

+d ( 4 ) * (uc ( sr-1 , sc-1 ) - uc ( sr-1 , 2 ) ) . . . 

+d(5) * (uc(sr,2) - uc ( sr , sc-1) ) . . . 

+ ( (2*dts A 2) /dxs) * (2*epss*time(l) + epss*all ( 1 , 1) ) ; 



vn(sr,sc) = vn(sr,l); 

% the compatability condition. 

innpd = (- (2*dxs) / (dts A 2) ) * (vn (sr , : ) - 2*vc(sr,:) +vo(sr,:)) 



+ bll ( 1 , :) ; 
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function cwv4 - coupled waves version 4 
written by: Lt. Hugh Me Bride USNR 
date : Apr 92 

cwv4 is the main program which couples the behaviour of the 'I-beam' shaped 
solid medium with the fluid medium. The elastic wave equation and the boundary 
conditions ,periodic and traction free are are satisfied for the solid, while 
the scalar wave equation and the periodic and radiating boundary conditions ar 
applied to the fluid. 



The general steps of the program are as follows 

1. The basic parameters are determined , that is the size of the domains 

as passed to it by the driver program cmodck.m 

2. All the global variables are calculated for both the fluid and the solid 
including the matrices required for the radiation boundary condition. 

3. vl and v3 are vectors used to find the amplitude, vl is the 
weighting vector 1/2 111.... 1 1/2 for the trapezoidal rule 
and v3 = exp ( i*n*pi*x) , the othhogonal vector. 

4. The I beam is built by basel 

5. The initial values of the u and v for the solid (un, vn , uc, vc , uo, vo) 

and m (mn me and mold) for the fluid are set to zero. 

6. The elastic wave equation and the periodic boundary conditions 
are applied to u and v 

7. The boundary conditions for the solid are applied 

8. The fluid and solid are coupled 

9. The freespace portions of the I beam are zeroed out 

10. The scalar wave equation and the periodic boundary conditions for 
the fluid are applied. 

11. The values of the amplitude are calculated and accumulated 

12. u, v and m are updated , that is the new values become 

the current values and the current values become the old values. 



function [v,un, vn,mn, z2 ] = cwv4 (m, kl , kt , kf , dts , dxs , x , nn, st , epss , tt , nt , mm) ; 
% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcoll and pclor are both required as there are two regions, one on 
% either side of the center spar of the I -beam which require 
% our attention and they contain the location of the nodes in question 
% kl - scaled longitudinal speed 
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% kt - scaled transverse speed 
% dxs - scaled spacing 
% dts - scaled time step 
% x - vector of length of the interval 
% nn - no of timesteps 

% st - stopping criteria for radiation boundary condition 
% epss - ratio of solid density to fluid density 
% nt - no of intervals in domain 

% mm - prpoagating mode under investigation ,mm = 0 fundamental 
% mm = 1 first and so forth 
% 2 ? - amplitude of propagating mode 



% Step 1 : r,c the dimensions of the domain, r-1 , c-1 and n-1 
% the dimensions of the interior 

[r c] = size (m) ; 
n = r-1; r = r - 1 ; c = c-1 ; 
axis ( 'xy ' ) 

% Step 2: The global variables 



% For the solid 

rho = rhov (kl , kt , dxs , dts) ; 
cunx = ucnx (kl , kt , dxs , dts) ; 
cuny = ucny (kl , kt , dxs , dts) ; 
cupx - ucpx (kl , kt , dxs , dts) ; 
cupy = ucpy (kl ,kt,dxs,dts) ; 

% For the fluid 

rhof =(dts~2) / (kf ~2*dxs~2) ; 
rcl = (2*dts~2) / (kf ~2*dxs A 2) 



ccs = (2 - 2*rho ( 1 ) -2 *rho ( 2 ) ) ; 
cvnx = vcnx (kl , kt , dxs , dts) ; 
cvny = vcny (kl , kt , dxs , dts ) ; 
cvpx = vcpx(kl ,kt,dxs,dts) ; 
cvpy = vcpy (kl ,kt,dxs,dts) ; 



ccf = (2 - 4*rhof) ; 



% The matrices for the radiation boundary condition allowing the propagating 
% modes to pass through the artificial boundary. 



anew = zeros(c+l); 
for pm = 0:st 

acurr = rbc (r+i , c+1 , pm, dxs , kf ) ; 
anew = anew+acurr; 

end 

ml = ( eye ( c+1 ) + (dts/ ( 2 *kf ~ 2 ) ) *anew) ; ml = inv(ml) ; 
m2 - (eye(c+l) - (dts/ ( 2 *kf * 2 ) ) *anew) ; 
t = trid (dxs , dts , c+1 , kf ) ; 

% Step 3: vl and v3 calculated so as to be able to determine the amplitude 
vl = del2 ( 2*nt ) ; 
v3 = exp ( i*mm*pi*x' ) ; 



% Step 4: basel builds the I -beam pclor,pcoll , coll, coir and fill 
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% allow us to identify the boundaries of the solid domain and the 
% corner nodes. 

[b,rows, pcoll,pcolr, fill ,coll , coir] = basel(n/2); 



% Step 5: All values for the fluid and the solid are initially 
% set to zero. 



un = 


zeros (size (ra) ) ; 


uc = un; 


uo = un; 


vn = 


zeros (size (m) ) ; 


vc = vn; 


vo = vn; 


mn = 


zeros (size (m) ) ; 


me = mn ; 


mold = mn 



for k - l:nn 

% Step 6: The elastic wave equation is applied to the interior 
% of the solid 

un(2:n,2:n) = ccs*uc (2 : n, 2 : n) - uo(2 :n, 2 :n) • . . 

+ rho(2) * (uc(l:n-l,2:n) + uc ( 3 : n+1 , 2 : n) ) . . . 

+ rho(l) * (uc(2 :n, l:n-l) + uc (2 : n, 3 : n+1) ) . . . 

+ rho(3 ) * (vc (1 :n-l , 1 : n-1) + vc (3 : n+1 , 3 :n+l) ) . . . 

- rho(3) * (vc(l:n-l, 3 :n+l) + vc ( 3 : n+1 , 1 : n-1 ) ) ; 



vn(2:n,2:n) = ccs*vc (2 : n, 2 : n) - vo(2 : n, 2 : n) . . . 

+ rho ( l) * (vc ( 1 : n-1 , 2 : n) + vc (3 : n+1 , 2 : n) ) . . . 

+ rho ( 2 ) * ( vc ( 2 : n, 1 : n-1) + vc ( 2 : n , 3 : n+1) ) . . . 

+ rho ( 3 ) * (uc ( 1 : n-1 , 1 : n-1) + uc ( 3 : n+1 , 3 : n+1) ) . . . 
- rho(3) * (uc(l:n-l, 3 :n+l) + uc ( 3 : n+1 , 1 : n-1) ) ; 



% The periodic boundary conditions for the solid 

un(2:n,l) = ccs*uc(2:n,l) - uo(2:n,l)... 

+ rho (2) * (uc ( 1 : n-1 , 1) + uc ( 3 : n+1 , 1) ) . . . 

+ rho(l) * (uc(2 :n, 2) + uc(2:n,n) ) . . . 

+ rho ( 3) * (vc ( 1 : n-1 , n) + vc ( 3 : n+1 , 2 ) ) . . . 

- rho ( 3 ) * ( vc ( 1 : n-1 , 2 ) + vc(3:n+l,n)) ; 

un(2:n,n+l) = un(2:n,l); 



vn(2:n,l) = ccs*vc (2 : n, 1) - vo(2:n,l)... 

+ rho ( 1) * (vc ( 1 : n-1 , 1) + vc (3 : n+1 , 1) ) . . . 

+ rho(2) * (vc(2 :n, 2) + vc (2 : n, n) ) . . . 

+ rho (3 ) * (uc ( 1 : n-1 , n) + uc (3 : n+1 , 2) ) . . . 
- rho ( 3 ) * (uc ( 1 : n-1 , 2) + uc(3:n+l,n) ) ; 

vn ( 2 : n # n+1) = vn ( 2 : n , 1 ) ; 



% Step 7: The boundaries of the I beam and including the corner nodes 
% are treated. 



102 



[un ,vn] = bnd f nx ( uc , vc , uo , vo , un , vn , rows , pcol 1 , cunx , cvnx ) ; 

[un ,vn] = bndfny (uc, vc,uo, vo,un, vn, rows, pcoll , pcolr , cuny , cvny) 
[un ,vn] = lcorny ( uc ; vc ; uo , vo , un , vn , rows , pcol 1 , cuny , cvny ) ; 

[un ,vn] = bndfpx(uc, vc,uo, vo , un , vn , rows , pcolr , cupx , cvpx ) ; 

[un ,vnj = bndfpy (uc, vc, uo,vo,un,vn, rows, pcoll, pcolr ,cupy,cvpy) 
[un ,vn] = lcorpy (uc, vc, uo,vo,un, vn, rows, pcoll, cupy,cvpy) ; 
[un,vn] = tcorl3 (uc, vc,uo, vo,un, vn, rows , pcoll , pcolr ,rho) ; 
[un,vn] = tcor24 (uc, vc,uo, vo,un, vn, rows , pcoll , pcolr ,rho) ; 



% Step 8 : The solid and the fluid are coupled (Note the is where the forcing 
% function of the problem is contained ) 

all = mc(l,:); bll = mc(2,:); 

[un ,vn,mnpd3 = coupl (uc , vc , uo , vo , un, vn , cupy , cvpy , k, dts , dxs , all , bll , epss) ; 



% Spep 9: Both sides of the center spar for all values of u and v are zeroed 
% out so as to prevent pollution 



un (rows , coll ) 
uc (rows, coll) 
un (rows, coir) 
uc(rows, coir) 
uo( rows, coll) 
uo (rows , coir) 



fill; vn (rows, coll) 
fill; vc (rows, coll) 
fill; vn(rows, coir) 
fill; vc (rows, coir ) 
fill; vo (rows , coll) 
fill; vo (rows , coir ) 



fill; 

ful- 

fill; 

fill; 

fill; 

fill; 



% Step 10: The scalar wave equation for the fluid 



% First the fluid/solid interface 

mn(l,2:c) = rhof * (mnpd ( 2 : c) + me (2 , 2 : c) . . . 
+mc(l,l:c-l) + me ( 1, 3 : c+1) ) +ccf *mc ( 1, 2 : c) . . . 
-mold ( 1 , 2 : c) ; 

% and it's periodic boundary condition 

mn(l,l) = rhof * (mnpd ( 1 , 1 ) + mc(2,l)... 
+mc(l,c) + mc(l,3)) +ccf *mc ( 1 , 1) . . . 

-mold (1,1) ; 

mn(l,c+l) = mn(l,l); 



% The interior points of the fluid 

mn(2:r,2:c) = rhof * (me ( 1 : r-1 , 2 : c) + me ( 3 : r+1 , 2 : c) . . . 
+mc ( 2 : r , 1 : c-1 ) + me ( 2 : r , 3 : c+1 ) ) +ccf *mc (2 : r , 2 : c) . . . 
-mold (2 : r , 2 : c) ; 



% The radiation boundary condition 

% Note: We are required to multiply a matrix by a row vector, but to do this 
% it must be transposed to a column . Matlab takes the Hermitian transpose 
% by default , so to ensure the correct signs we must void this effect by 
% taking the conjugate of the transpose before we do our calculations. 

% The process is then reversed so as the updated value has the correct dimension 
% 
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si = conj (me (r, : ) ' ) ; 
s2 = con j (me (r+1 , : ) ' ) ; 
s3 = conj (mold (r+1 ; 
int = ml*(rcl*sl + t*s2 - m2*s3) ; 
mn(r+l,:) = conj(int'); 



% the periodic boundary condition for the artificial boundary 

mn(2:r,l) = rhof * (mc(l :r-l, 1) + mc(2:r,2) + ... 

me (3: r+1 ,1)+ mc(2:r,c)) + ccf *mc (2 : r , 1) - mold(2:r,l); 

mn(2:r,c+l) = mn(2:r / l); 



% Step 11: Calulates amd accumulates the values of the amplitude 
aa = abs ( (mn (r , : ) . *vl ) *v3 ) ; 

22 = [ z2 aa]; 

% Step 12: The values for u,v and m are updated. 

uo = uc; uc = un; 

vo = vc; vc = vn; 

mold = me; me — mn; 



end 
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% function del2 - del funtion used for thr trapezoidal rule 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% bulids a vector of appropiate length used to weight the elements of the 
% quantity being integrated 

function d = del2(n) 
a = ones ( 1 , n+1) ; 
a(l) = .5; a (n+1) = .5; 
a2 = ones ( 1 , n+1) /n; 

d = a.*a2; 
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% function lcorny - left corner facing in the negative y direction 
% i.e unit normal = (0,-1) 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% 

% 

% lcorny performs the same function as bndfny except it only operates on 
% 4 points , those which lie at the extremities of the domain, the points 
% as if fig( ) Periodicity is used for pseudonodes which lie outside the domain. 
% This requires us to pick out each nodes individually (9 for each 4 points , 

% making a total of 36 ) they are weighted in the same fashion 

% as in bndfny and the updated values for u and v are calculated for both sides 

function [un ,vn] = lcorny (uc , vc , uo , vo , un , vn , rows , pcoll , c , d) ; 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcoll allows us to pick out the column location of those nodes 
% to the left of centerline and using the periodicity we substitute 
% this value into the corresponding point on the right of 
% centerline. 



% c = cuny; 
% d = cvny; 



[ i 1 , j 1 ] = size(pcoll) ; 
[ i 2 , j 2 ] = size(rows); 

[ i3 , j 3 ] = size (uc) ; 



% the updated values of u to the left of centerline are 
% calculated 

ucc = uc (rows (j2 ) +1 , pcoll ( 1) ) ; 
uccl = uc (rows (j2) +2 , pcoll (1) ) ; 
ucr = uc (rows ( j2 ) +1 , pcoll (2 )) ; 
ucl = uc (rows ( j2) +1 , j3-l) ; 

ver = vc (rows (j2 ) +1 , pcoll (2 )) ; 
vru = vc ( rows ( j 2 ) +2 , pco 11(2)); 
vlu = vc (rows ( j2 ) +2 , j3-l) ; 
vcl = vc (rows ( j2 ) +1 , j 3 — 1 ) ; 

uo 1 = uo ( r ows ( j 2 ) + 1 , pco 1 1 ( 1 ) ) ; 

ul = c(l)*ucc - uol + c(2)*uccl +c(3)*(ucr + ucl)... 
+c(4)*(vlu - vru) + c(5) * (ver- vcl); 



uccl = uc (1 , pcoll (1) ) ; 
ucul = uc (2 , pcoll (1) ) ; 
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ucrl = uc(l,pcoll ( 2 ) ) ; 
ucll = uc ( l , j 3-1) ; 

vcrl = vc (1 ,pcoll ( 2) ) ; 
vrul = vc (2 , pcoll ( 2 ) ) ; 
vlul = vc ( 2 , j 3-1) ; 
veil = vc(l, j3-l) ; 

uoll = uo( 1, pcoll (1) ) ; 



ull = c(l)*uccl - uoll + c(2)*ucul +c(3) * (ucrl + ucll)... 
+c(4)*(vlul - vrul) + c(5)*(vcrl- veil); 



% the updated values of v to the left of centerline are 
% calculated 



vucc = vc (rows ( j2) +1 , pcoll ( 1) ) ; 
vuccl = vc(rows ( j2) +2 , pcoll (1) ) ; 
vucr = vc(rows ( j2) +l,pcoll (2) ) ; 
vucl = vc (rows ( j2) +1, j3-l) ; 

wer = uc (rows ( j 2 ) +1 , pcoll ( 2 ) ) ; 
wru = uc (rows ( j 2 ) +2 , pcoll ( 2 ) ) ; 
wlu = uc (rows ( j2) +2 , j 3 — 1 ) ; 
wcl = uc(rows ( j2) +1, j3-l) ; 

vuol = vo ( rows ( j 2 ) +1 , pool 1(1)) ; 

vul = d(l)*vucc - vuol + d(2)*vuccl + d(3)*(vucr + vucl)... 
+d (4 ) * (wlu - wru) + d(5)*(vvcr- wcl) ; 



vuccl = vc ( 1, pcoll (1) ) ; 
vucul = vc (2 , pcoll ( 1) ) ; 
vucrl = vc ( 1 , pcoll (2) ) ; 
vucll = vc(l,j3-l); 

wcrl = uc ( 1, pcoll (2) ) ; 
vvrul = uc (2 , pcoll ( 2) ) ; 
wlul = uc(2, j 3 - 1 ) ; 
well = uc(l,j3-l); 

vuoll = vo ( 1 , pcoll (1)); 



vull = d(l) *vuccl - vuoll + d(2)*vucul +d(3)*(vucrl + vucll)... 
+d(4) *(wlul - vvrul) + d(5)*(wcrl- well); 

% the values put in the correct positions 

% and since we have periodic boundary conditions, we insert the left 
% hand value into the correspondding right hand position 

un (rows ( j2) +1 , pcoll ( 1) ) = ul; 
un(rows( j2) +1 , j3) = ul; 
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vn (rows (j2) +1, pcoll (1) ) = vul; 
vn (rows ( j2) +1, j3) = vul; 



un (1, pcoll (1) ) = ull; 
un(l,j3) = ull; 

vn (l, pcoll ( 1) ) = vull; 
vn(l,j3) = vull; 
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% function lcorpy - left corner facing in the positive y direction 
% i.e unit normal = (0,1) 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% 

% 

% lcorpy performs the same function as bndfpy except it only operates on 
% 2 points , those which lie at the extremities of the domain, the points 
% as in fig( ) Periodicity is used for pseudonodes which lie outside the domain. 
% This requires us to pick out each nodes individually (9 for each 2 points , 

% making a total of 18 ) they are weighted in the same fashion 

% as in bndfpy and the updated values for u and v are calculated for both sides 



function [un ,vn] = lcorpy ( uc , vc , uo , vo , un , vn , rows , pcoll , c , d) ; 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rows : used to identify the elements of the matrix which 

% are zero for all times they also contain the row location of 
% the nodes on the boundary. 

% pcoll allows us to pick out the column location of those nodes 
% to the left of centerline and using the periodicity we substitute 
% this value into the corresponding point on the right of 
% centerline. 



% c = cupy; 

% d = cvpy; 



% there are now only two values which need to be calculated 
% as the other boundary facing in the positive 
% y direction is at the fluid/solid interface and requires 
% a special treatment 

[il,jl] = size(pcoll); 

[i2,j2] = size (rows); 

[ i3 , j 3 ] = size (uc) ; 



% the updated u value 

ucc = uc (rows ( 1) -1 , pcoll (1) ) ; 
ucu = uc (rows ( 1 ) -2 , pcoll ( 1 )) ; 
ucr = uc (rows ( 1) -I , pcoll (2 )) ; 
uc 1 = uc (rows ( 1) -1 , j3-l) ; 

ver = vc (rows ( 1) -1 , pcoll (2) ) ; 
vru = vc (rows ( 1 ) -2 , pcoll (2 ) ) ; 
vlu = vc (rows ( 1) -2 , j3-l) ; 
vcl = vc(rows(l) -1, j3-l) ; 

uol = uo (rows ( 1 ) -1 , pcoll ( 1) ) ; 

ul = c(l)*ucc - uol + c(2) *ucu +c(3)*(ucr + ucl) . . . 
+c ( 4 ) * (vlu - vru) + c (5) * (ver- vcl); 
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% the updated v value 



vucc = vc (rows (1) -1 , pcoll ( 1) ) ; 
vucu = vc (rows ( 1) -2 , pcoll ( 1 ) ) ; 
vucr = vc (rows (1) -1 , pcoll (2) ) ; 
vucl « vc (rows (1) -1 , j3-l) ; 

wcr = uc (rows (1) -1 , pcoll (2 )) ; 
wru = uc (rows ( 1) -2 , pcoll (2 ) ) ; 
wlu = uc (rows ( 1) -2 , j3-l) ; 
wcl = uc(rows (1) -1, j 3 — 1 ) ; 

vuol = vo (rows (1) -1 , pcoll ( 1) ) ; 

vul = d(l)*vucc - vuol -i- d(2)*vucu + d(3)*(vucr + vucl)... 
+d(4)*(wlu - wru) + d(5)*(wcr- wcl) ; 



% using periodicity we update the values to the left and right of 
% the center spar 

un (rows ( 1) -1 , pcoll (1) ) = ul; 
un (rows ( 1) -1 , j 3 ) = ul; 
vn (rows ( 1) -1 , pcoll (1) ) = vul; 
vn ( rows ( 1 ) - 1 , j 3 ) = vul; 
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% function rbc - radiating bondary condition 
% written by: Lt. Hugh Me Bride USNR 
% date : Mar 92 

% 

% rbc builds the matrix A(i,k) which is required by the radiation 
% boundary condition 

function v = rbc ( j , 1 , pm, dxs , kf ) 

% variables 

% j,l -dimensions of the matrix 
% pm - propagating modes 
a = zeros ( j ) ; 

for i = 1 : j , 

for k = 1:1, 

a(i,k) = exp(sqrt (-1) *pm*pi* (i-k) *dxs) ; 
end 



end 



a ( : , 1) = • 5*a ( : , 1 ) ; a(:,k) = .5*a(:,k); 

v = sqrt (kf *2- (pm*pi) *2 ) *a ; 
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% function rho- del funtion used for thr trapezoidal rule 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% generates a vector containing constants used repeatedly throughout 
% the program. 

function rho = rhov (kl , kt , dxs , dts) ; 



rhol =(dts A 2) / (kl A 2*dxs A 2) ; 
rho2 = (dts A 2) / (kt A 2*dxs A 2) ; 

rho3 = (dts A 2/ ( 4 *dxs A 2 ) ) *( (l/kl A 2)-(l/kt A 2) ) ; 
rho = [rhol rho2 rho3]; 
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% function stop - defines the stopping criteria for constructing the radiation 
% boundary condition i.e if the result of stop is 0 then only the fundamental mo 
% is allowed to propagate, 1 only the fundamental and the first mode are allowed 
% to prpogate and so on. 

% written by: Lt. Hugh Me Bride USNR 
% date : Mar 92 

% 

function v = stop(k) 



% variables 
% k - 

% n - nth propagating mode 
n = 0 ; m = -1; 
while (k A 2 - n~2*pi A 2) > 0. 
m = m+1; n = n+1; 

end 
v - m; 
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% function tcorl3 - treatment of corners 1 and 3 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% 

% tcor!3 applies the elastic wave equation at 
% corners 1 and 3 as per fig( ) since they use the 
% same difference formula. 

% the corner nodes are located (1 first then 3) 

% identified as p and q (pi , ql in the case of corner 3) 

% the necessary neighbours are picked off from the u and v 
% matrices and weighted accordingly and the new updated values 
% for u and v are computed. 



function [un,vn] = tcor 13 (uc, vc, uo, vo , un , vn, rows , pcoll , pcolr , rho) 



% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rho : vector containing global constants 

% rows : used to identify the elements of the matrix which 
% are zero for all times they also contain the row location of 
% the corner nodes. 

% pclor and pcoll carries out the same function as rows for 
% the columns, and the contain the column location of the 
% corner nodes. 



[il,jl] = size(rows); 

[ i2 , j 2 ] = size(pcolr); 

( i3 , j 3 ] *= size (pcoll) ; 

% we pick off the elments of rows and pcolr which 
% identify the location of corner 1. 

p = rows(il)-l; 
q = pool 1 ( j 3 ) ; 

% generate any required local constants 

cc = (2 - 2*rho(l) -2*rho(2) ) ; 

% the updated values of the corner nodes for u and v are 
% calculated 



un ( p, q) = cc*uc (p , q) - uo(p,q)... 

+ rho(l) *(uc(p,q-l) + uc (p , q+1) ) . . . 

+ rho(2) * (uc(p-l,q) + uc(p+l,q) ) . . . 

- 2*rho(3) * (vc(p-l ,q-l) + vc(p+l,q+l) + 2 *vc (p , q) ) . . . 

+ 2*rho(3) *(vc(p,q+l) + vc(p,q-l) + vc(p+l,q) + vc(p-l,q)); 



vn(p,q) = cc*vc(p,q) - vo(p,q)... 

+ rho ( 2 ) * ( vc (p, q-1) + vc (p, q+1) ) . . . 

+ rho ( 1 ) * ( vc ( p-1 , q) + vc (p+1 , q) ) . . . 
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- 2*rho (3 ) * (uc (p-1 ,q-l) + uc(p+l,q+l) + 2*uc(p,q) ) . . . 

+ 2*rho(3) * (uc(p,q+l) + uc(p,q-l) + uc(p+l,q) + uc(p-l,q)); 



% the same process is repeated for corner 4 below 



pi = rows(jl)+l; 
ql = pcolr (i2) ; 



un(pl,ql) =cc*uc (pi , ql) - uo(pl,ql)... 

+ rho(l) * (uc(pl ,ql-l) + uc (pi , ql+1) ) . . . 

+ rho(2) * (uc (pl-1 ,ql) + uc (pl+1 , ql ) ) . . . 

- 2 *rho (3) *(vc(pl-l,ql-l) + vc (pl+1 , ql+1 ) + 2*vc (pi , ql) ) . . . 

+ 2*rho(3) * (vc (pi , ql+1) + vc(pl,ql-l) + vc(pl+l,ql) + vc (pl-1 ,ql) ) ; 



vn(pl,ql) =cc*vc (pi ,ql) - vo(pl,ql)... 

+ rho (2) * (vc (pi, ql-1) + vc (pi , ql+1 ) ) . . . 

+ rho (1) * (vc (pl-1 ,ql) + vc (pl+1 , ql) ) . . . 

- 2*rho (3 ) * (uc (pl-1 ,ql-l) +uc(pl+l ,ql+l) +2*uc (pi , ql) ) . . . 

+ 2*rho ( 3 ) * (uc (pi , ql+1) + uc(pl,ql-l) + uc(pl+l,ql) + uc (pl-1 , ql ) ) ; 
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% function tcor24 - treatment of corners 2 and 4 
% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% 

% 

% tcor24 applies the elastic wave equation at 
% corners 2 and 4 as per fig( ) since they use the 
% same difference formula. 

% the corner nodes are located (2 first then 4) 

% identified as p and q (pi ,ql in the case of corner 4) 

% the necessary neighbours are picked off from the u and v 
% matrices and weighted accordingly and the new updated values 
% for u and v are computed. 



function [un,vn] = tcor24 (uc , vc , uo , vo , un , vn , rows , pcoll , pcolr , rho) ; 
% variables 

% un : updated values of u vn : updated values of v 

% uc : current values of u vc : current values of v 

% uo : old values of u vo : old values of v 

% rho : vector containing global constants 

% rows : used to identify the elements of the matrix which 
% are zero for all times they also contain the row location of 
% the corner nodes. 

% pclor and pcoll carries out the same function as rows for 
% the columns, and the contain the column location of the 
% corner nodes. 



[il,jl] = size (rows) ; 

[i2,j2] = size(pcolr); 

[i3,j3] = size(pcoll); 

% we pick off the elments of rows and pcolr which 
% identify the location of corner 2. 

p = rows ( il ) -1 ; 
q = pcolr ( i 2 ) ; 

% generate any required local constants 
cc = (2 - 2 *rho ( 1) -2 *rho ( 2 ) ) ; 



% the updated values of the corner nodes for u and v are 
% calculated 

un (p,q) = cc*uc(p,q) - uo(p,q)... 

+ rho(l) * (uc(p,q-l) + uc (p, q+1 ) ) . . . 

+ rho (2) * ( uc (p-1 , q) + uc (p+1 , q) ) . . . 

+ 2*rho(3) * (vc(p+l,q-l) + vc(p-l,q+l) + 2 *vc (p, q) ) . . . 

- 2*rho(3) * (vc(p,q+l) + vc(p,q-l) + vc(p+l,q) + vc(p-l,q)); 



vn (p , q) = cc*vc(p,q) - vo(p,q)... 

+ rho(2) * (vc (p, q-1) + vc(p,q+l) ) . . . 

+ rho(l) * (vc (p-1 , q) + vc (p+ 1 , q) ) . . . 

+ 2*rho(3) * (uc(p+l,q-l) + uc(p-l,q+l) + 2*uc (p, q) ) . . . 

- 2*rho(3) *(uc(p,q+l) + uc(p,q-l) + uc(p+l,q) + uc(p-l,q)); 
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% the same process is repeated for corner 4 below 



pi = rows(jl)+l; 
ql= pcol 1 ( j 3 ) ; 



un (pi , ql ) = cc*uc(pl,ql) - uo(pl,ql)... 

+ rho ( 1) * (uc (pi # ql-l) + uc (pi , ql+1 ) ) . . . 

+ rho (2) *(uc(pl-l,ql) + uc (pl+1 , ql) ) . . . 

+ 2*rho(3) * (vc (pl+1 , ql-1) + vc (pl-1 , ql+1) + 2*vc (pi , ql) ) . . . 

- 2*rho ( 3 ) * (vc (pi , ql+1) + vc(pl,ql-l) + vc(pl+l,ql) + vc(pl-l,ql)) 



vn (pi , ql) = cc*vc(pl,ql) - vo(pl,ql)... 

+ rho (2) *(vc(pl / ql-l) + vc (pi , ql+1 ) ) . . . 

+ rho(l) * (vc(pl-l,ql) + vc (pl+1 , ql) ) . . . 

+ 2* rho (3)*(uc( pl+1, ql-1) + uc (pl-1 , ql+1) + 2 *uc (pi , ql ) ) . . . 

- 2*rho (3) *(uc(pl,ql+l) + uc(pl,ql-l) + uc(pl+l,ql) + uc (pl-1 , ql ) ) ; 
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% function trid - tridiagonal 
% written by: Lt. Hugh Me Bride USNR 

% date : Mar 92 

% 

% trid generates a tridiagonal matrix for the 
% radiation boundary condition of the fluid . 

% the elements of the sub and super diagonal are (dts/kf*dxs) *2 
% the main diagonal component is 2 - 4 (dts/kf *dxs) "2 
% the (n # 2) and (l,n-l) contain (dts/kf *dxs) A 2 to satisfy the 
% periodic boundary conditions. 

function m = trid (dxs, dts , n , kf ) 

% variables 

% dxs : scaled spacing 

% dts : scaled time step 

% n : dimension of matrix 

% kf : scaled constant 



% an identity matrix for the elements of the main diagonal 
dl = eye (n) ; 

% the sub and super diagonals 

d2 = diag (ones (n-1 , 1) , 1) + diag (ones (n-1, 1) , -1) ; 
% the required co-efficients 

rho = dts/ (kf *dxs) ; rho = (rho~2) ; 

% generates the required matrix 

d =2*dl -4*rho*dl +rho*d2; 
d(l,n-l)= rho ; d(n,2) = rho; 



m 



= d; 
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% function ucnx - u coefficients/constants for the boundary facing the 
% negative x direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 

% technique to the boundary with unit normal (-1 , 0) for the corresponding u 

% and v values. 

% 

function cunx = ucnx (kl , kt , dxs , dts) 

% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 

% dts - scaled time step 

% the coefficients are calculated and stored in the vector cunx 
% for use in bndfnx 



cl = 2 - 2 * ( dts A 2 /dxs A 2 ) * ( (l/kl*2) + ( 1/kt A 2 ) ) ; 
c2 = 2*dts A 2/ (dxs~2*kl~2) ; 
c3 = (dts~2) / (dxs A 2*kt~2) ; 

c4 = - (dts A 2/ (2*dxs A 2) )*( (l/kl~2) -(l/kt A 2) ) ; 
c5 = (dts^2/ (2 *dxs A 2 ))*( (l/kl A 2)-(3/kt A 2) ) ; 
cunx = [cl c2 c3 c4 c5] ; 
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% function ucny - u coefficients/constants for the boundary facing the 
% negative y direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 
% technique to the boundary with unit normal (0 , -1 ) for the corresponding u 
% and v values. 

% 

function cuny = ucny (kl , kt , dxs , dts) 



% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 
% dts - scaled time step 

% the coefficients are calculated and stored in the vector cuny 
% for use in bndfny 

cl = 2 - 2* (dts A 2/dxs A 2 ) * ( (l/kl A 2) +(l/kt~2) ) ; 
c2 = 2* (dts A 2/ (dxs A 2*kt A 2) ) ; 
c3 = (dts A 2 ) / (dxs A 2 *kl A 2 ) ; 

c4 = - (dts A 2 / (2*dxs*2) ) * ( ( 1/kl ~2 ) - ( l/kt A 2 ) ) ; 
c5 = ( dts A 2 / ( 2 *dxs A 2 ) ) * ( ( 3 /kt A 2 )-(l/kl A 2)) ; 
cuny = [cl c2 c3 c4 c5] ; 
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% function ucpx - u coefficients/constants for the boundary facing the 
% positive x direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 

% technique to the boundary with unit normal (1,0) for the corresponding u 

% and v values. 

% 

function cupx = ucpx (kl , kt , dxs , dts) 

% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 

% dts - scaled time step 

% the coefficients are calculated and stored in the vector cupx 
% for use in bndfpx 



cl = 2 - 2 * (dts * 2 /dxs A 2 )*((l/kl A 2)+( 1/kt ~2 ) ) ; 
c2 = 2* (dts A 2/ (dxs A 2*kl*2) ) ; 
c3 = (dts~2) / (dxs~2*kt~2) ; 
c4 = (dts A 2/ (2*dxs~2) ) *( ( 1/kl A 2 ) - ( 1/kt " 2 ) ) ; 
c5 = - (dts*2/ (2*dxs A 2) ) * ( (l/kl*2) -(3/kt"2) ) ; 
cupx = [cl c2 c3 c4 c5] ; 
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% function ucny - u coefficients/constants for the boundary facing the 
% positive y direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 

% technique to the boundary with unit normal (0,1) for the corresponding u 

% and v values. 

% 

function cupy = uepy (kl , kt , dxs , dts) 

% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 

% dts - scaled time step 

% the coefficients are calculated and stored in the vector cupy 
% for use in bndfpy 



cl = 2 - 2* (dts~2/dxs A 2) * ( ( 1/kl A 2 ) + ( 1/kt A 2 ) ) ; 
c2 = 2* (dts~2/ (dxs A 2*kt A 2 ) ) ; 
c3 = (dts A 2 ) / (dxs A 2 *kl A 2 ) ; 

c4 = (dts A 2 / ( 2*dxs*2 ) )*((l/kl~2)-( l/kt*2 ) ) ; 
c5 = (dts A 2/ (2*dxs A 2) ) * ( ( l/kl A 2) - (3/kt A 2) ) ; 
cupy = [cl c2 c3 c4 c5] ; 
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% function vcnx - v coefficients/constants for the boundary facing the 
% negative x direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 
% technique to the boundary with unit normal (-1 , 0) for the corresponding u 
% and v values. 

% 



function cvnx = vcnx (kl , kt , dxs , dts) 



% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 
% dts - scaled time step 

% the coefficients are calculated and stored in the vector cvnx 
% for use in bndfnx 



dl = 2 - 2 * (dts A 2 / dxs A 2 ) *( (l/kl A 2) + ( 1 /kt A 2 ) ) ; 
d2 = 2* (dts A 2/ (dxs A 2*kt A 2) ) ; 
d3 = (dts A 2) / (dxs A 2 *kl * 2 ) ; 

d4 = -(dts"2/ (2*dxs*2) ) * ( (l/kl A 2) -(l/kt~2) ) ; 
d5 = ( dts A 2/ ( 2 *dxs A 2 ) ) * ( ( 3 /kt * 2 )-(l/kl*2)) ; 
cvnx = [dl d2 d3 d4 d5]; 
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% function vcny - v coefficients/constants for the boundary facing the 
% negative y direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 
% technique to the boundary with unit normal ( 0 , -1 ) for the corresponding u 
% and v values. 

% 



function cvny = vcny (kl , kt , dxs , dts) 

% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 
% dts - scaled time step 

% the coefficients are calculated and stored in the vector cvny 
% for use in bndfny 



dl = 2 - 2*(dts^2/dxs^2) *( ( 1/kl A 2 ) + ( l/kt^2 ) ) ; 
d2 = 2* (dts^2/ (dxs^2*kl^2) ) ; 
d3 = (dts~2) / (dxs~2*kt~2) ; 

d4 = -(dts*2/ (2*dxs~2) )*( ( 1/kl *2 ) - ( l/kt*2 ) ) ; 
d5 = (dts^2 / (2 *dxs*2 ) ) * ( ( l/kl*2 ) - ( 3 /kt A 2 ) ) ; 
cvny = [dl d2 d3 d4 d5]; 



124 



% function vcpx - v coefficients/constants for the boundary facing the 
% positive x direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 

% technique to the boundary with unit normal (1,0) for the corresponding u 

% and v values. 

% 

function cvpx = vcpx (kl , kt , dxs , dts) 

% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 

% dts - scaled time step 

% the coefficients are calculated and stored in the vector cvpx 
% for use in bndfpx 



dl = 2 -2* (dts~2/dxs~2) * ( (l/kl / '2) + (l/kt^2) ) ; 
d2 = 2* (dts~2 / (dxs *2 *kt~2 ) ) ; 
d3 = (dts*2) / (dxs*2*kl~2) ; 

d4 = ( dts~ 2 / ( 2 *dxs A 2 ) )*( (l/kl A 2)-( 1/kt ~ 2 ) ) ; 
d5 = (dts'' 2/(2 *dxs A 2 ) )*( (l/kl~2)-(3 /kt~2 ) ) ; 
cvpx = [dl d2 d3 d4 d5]; 
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% function vcny - v coefficients/constants for the boundary facing the 
% negative y direction. 

% written by: Lt. Hugh Me Bride USNR 
% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 
% technique to the boundary with unit normal (0,1) for the corresponding u 
% and v values. 

% 



function cvpy = vepy (kl , kt , dxs , dts) 

% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 
% dxs - scaled spacing 
% dts - scaled time step 

% the coefficients are calculated and stored in the vector cvpy 
% for use in bndfpy 



dl = 2 - 2*(dts A 2/dxs~2)*( (l/kl A 2)+(l/kt A 2) ) ; 
d2 = 2* (dts A 2/ (dxs A 2*kl^2) ) ; 
d3 = (dts A 2) / (dxs*2*kt^2) ; 
d4 = (dts A 2 / ( 2 *dxs ^ 2))*((l/kl A 2)-( 1/kt A 2 ) ) ; 
d5 = - (dts A 2/ (2*dxs"2) ) * ( (l/kl"2) -(3/kt A 2) ) ; 
cvpy = [dl d2 d3 d4 d5]; 
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